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Polynomial Root Solving on the Electronic 
Differential Analyser (A Technique for Finding 
the Real and Complex Roots of a Polynomial 
using an Electronic Differential Analyser) 


1. Introduction. This paper presents a technique for the solution of high 
degree polynomials with real coefficients. By means of this technique the real 
and complex roots of the polynomial 


F(%) = aon + ayx™ + +++ +a? 


can be obtained. The technique is essentially a graphic one in which the function 
f(x) and its derivatives f’ (x), f’’ (x), --- are generated by an electronic differential 
analyser [1]. The method depends heavily on the ability of the analyser to pro- 
duce an accurate plot of the function and its derivatives. The use of these accurate 
graphs is obvious for finding the real roots of the equation. But a systematic use 
of graphs for the finding of complex roots is not so obvious and hence not widely 
used [2]. 

This computer approach to the finding of the complex roots requires addi- 
tional curves of the real and complex parts of the complex function 


U +4V = f(z) = aor* + aye” " +---+ an; (¢ = x + ty). 


The real part U(x, y) and the complex part V(x, y) are also plotted by the analyser 
and complex roots are found from the common zeros. 

2. Method of Solution. 

1. Real Roots. In order to generate the equation y = f(x) on the differential 
analyser the polynomial must be differentiated times, where n is the degree of 
the polynomial. This will produce the derivatives f*(x), f*-'(x), ---, f’(x) where 
f"(x) is a constant, say f*. If a voltage of value f* is introduced into an integrating 
amplifier the output is xf". (See Fig. 1a.) Putting the value of f*(0) into the 
same integrator makes the output f*'(x) (neglecting the sign change of the 
amplifier). This process is continued until the final integration produces f(x). 
The real roots are picked from the curve as the values of x for which f(x) = 0. 

2. Complex Roots. The technique for finding the complex roots of the poly- 
nomial is based on the idea of a surface f(z) over the horizontal xy plane where 


f(z) = f(x + ty) = U(x, y) + iV (x, y). 


For f(z) real, V(x, y) = 0; f(z) = U(x,¥) so U(x, y) is a space curve whose 
projection on the xy plane is V(x, y) = 0. 

For the solution of the polynomial equation f(z) = U + 2V = 0 both the real 
function U and the imaginary function V must be zero so the root x; + iy: 
represents the intersection of the space curve U(x, y) with the curve V(x, y) = 0. 
This is a point in the xy plane where the U function is zero. 
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Fic. 1. General wiring nay cag for polynomial root solving. 
Sign inverters have not been shown. 


The differential analyser can handle only one variable of integration so the 
common point in the xy plane where both U and V are zero must be approached 
by a judicious selection of one of the variables, usually y = y;. Then a plot of the 
two functions U(x, y;) and V(x, 1) reveals how close ¥,; is to the actual imaginary 
part of the root. If the plot of U; and V; indicates a common crossing of the axis, 
say at x;, then the chosen y, is the correct one and the root is x; + iy:. Notice 
that this also yields the conjugate root x; — 71. 

If the curves U; and V; do not have a common crossing of the abscissa, a 
change in y; is made and another plot of the two equations is generated. This is 
continued until the root is obtained to the desired accuracy. 

The judicious choice of y; is assisted by a knowledge of the positions of maxima, 
minima, points of inflection, and accessory points of inflection (where a maximum 
and a minimum have run together) of the real curve. 

A. J. KEMpNER’s [2] suggestion that a study of the derivatives is helpful in 
estimating the complex roots of a polynomial is borne out by the nature of the 
real and imaginary parts of the complex function f(z). The real part 


U(x, y) = f(x) — f’(x)y?/2! + FEY (x)yt/4! — «+, 

















POLYNOMIAL ROOT SOLVING ON ELECTRONIC DIFFERENTIAL ANALYSER 141 





f(x) 
x 
f 


‘ 





Fic. 2. Graph showing real roots of quintic equation x* + x* — 7x* — 22x* +x +1 =0of 
example 1. Note: Change of vertical scale at X = 1. 


and the imaginary part 


Vix, y) = f'(x)y — f'" (x)y8/3! + $M (x)y8/5! — ---; 


that is, the real part is the alternating sum of the polynomial and its even deriva- 
tives multiplied by the corresponding even powers of y and inverse factorials, 
while the imaginary part is the sum of the odd derivatives multiplied by the odd 
powers of y and corresponding inverse factorials. 

Since the respective derivatives are generated in order to produce the function 
f(x) (as shown in Fig. 1a) they are available for finding a first approximation to 
the imaginary part of a complex root. They are then used for producing U(x, ¥:) 
and V(x, ¥:) after being multiplied by the appropriate factor (e.g., f(x) is 
multiplied by y:?/2!). (See Fig. 1b.) 

3. Examples of Computer Solutions. 

Example 1. Real Roots of a Quintic Equation. The polynomial 


x5 + xt — 7x8 — 22x° + 24+1 = f(x) 


is used as an example for the machine solution of real roots. The roots sought are 
the real positive roots [3]. This requires only the generation of the respective 
derivatives (Fig. 1a) to produce f(x). The computer plotted the function as shown 
in Fig. 2. The curve crossed the x axis at the roots 0.22 and 3.28. 

Example 2. Complex Roots of a Quartic Equation. The polynomial 


xt + x? + Sx? + 5x + 12 = f(x) 


is used as an example for the machine solution of complex roots [4 ]. The computer 
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Fic. 3. Graph showing common intersection of x axis by U(x, 1) and V(x, y:) indicating the 
roots — 1 + iv2 of the equation x* + x* + 5x* + Sx + 12 = 0. 
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is wired as in Fig. 1a and Fig. 1b and U and V curves are plotted for negative x 
with y = y2. Fig. 3 indicated the simultaneous crossing of the axis at x = 1. 
Since the computer plotted the negative halves of the curves this crossing is then 
the conjugate pair — 1 + iy2. 

Example 3. Sixth Order Polynomial with Three Pairs of Complex Roots. The 
polynomial [5 ] 


x8 — x5 + 2xt — 3x? + 2x? +2441 = f(x) 


was chosen as a fairly difficult equation to solve by this method since all roots are 
complex and clustered about the origin. 

Figure 4 is a plot of the expressions U and V for the value of y; = 0.64. The 
pair of roots as given by the computer was 1.0 + 0.647. (These roots are given 
in USPENSKY as 1.08018 + 0.639037.) 
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Fic. 4. Graph showing common intersection of x axis , U(x, v1) and V(x, 1) nen Tyme the 
roots 1.08 + 0.64: of the equation x* — x® + 2x — 3x8 + 2x7 ++%+1=0. 


4. Conclusion. The method of root solving presented herein is a new technique 
employing the latest type of analog computer to the old problem of root solving. 
The solutions generated by the computer are fast as compared with the compu- 
tational methods or the mechanical root finders. The accuracy of roots can easily 
be kept in the range of engineering accuracy. Mathematical manipulations of the 
original equation are minimized by the method which requires only the derivatives 
of the polynomial whose roots must be determined. 


CyriL ATKINSON 
University of Calif. 
Berkeley, Calif. 
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Notations for Partitions 


The present paper gives a contribution to the automatic calculation of the 
characters of certain symmetric groups, as described in a paper which recently 
appeared in MTAC [1]. The repeated application of the recursion formula (7) of 
the cited paper involves the storing of a chain of partitions. The allotted storage 
space must be sufficient to meet even the case of a chain of as many as about 2/2 
partitions. In other words, provision must be made for storing about n?/4 integers, 
components of these partitions. It would be a waste of memory space to store 
these integers separated, numerous and small as they are. Furthermore, as 2 
increases, it will be still more important that the partitions should be available 
in the rapid access storage of the machine. Therefore, they should be stored by 
some method of packing, suitable for the operations to be applied, particularly 
during the calculation of characters. 

In the sequel, I shall briefly describe three such methods, which might perhaps 
be of some interest. The first method [2] was used in an experimental program 
for character calculation with the Swedish relay computer BARK (August, 1953) 
and also in the first version of the character program using the Swedish electronic 
computer BESK (Spring, 1954). The second method is employed in the actual 
version of this program (tested on BESK, June, 1954). The third method, a 
modification of the second, has not yet been used in the program, though it is 
described here for theoretical reasons. 

Preliminaries. Theoretically, the minimum of the number, x, of binary digits 
required to represent all the partitions of an integer m, is the smallest integer, not 
less than logs P(m), where P(m), the number of the partitions, can be estimated 
by means of Harpy-RAMANUJAN’S formula [3]. The result obtained, 


Xmin ~ TV2/3 log: e-Vn — logs 4nV3 ~ 3.70Vn — logs n — 2.79 


indicates an ideal towards which the packing methods should lead. 

If packing with fixed boundaries between the packed elements, one should 
store the partition components, each diminished by 1. In this case, m must also 
be stored. The value of x then lies between m and 3n/2. For example, in a 40 digit 
register one can store the partitions of m < 30 by this method. It will be seen that 
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our first and second methods yield slightly better values of x, namely x = n 
and x = n + 1, respectively, whereas our third method gives approximately 
x < V¥32n + 4 — 2, which is, for m > 28, in better accord with the Hardy- 
Ramanujan value. 

Our notations for partitions are closely related to the YouNG diagrams of the 
partitions [4], which are important in the NAKAYAMA version [5 ] of the recursive 
method. This version simplifies the calculation scheme, and it has, therefore, been 
employed in our program. 

It is essential, in the methods which we are going to describe, that the com- 
puting machine is equipped with the operation known as “normalization [6].” 
We assume that the machine represents numbers modulo 2 by WN digit words, 
the fractional part consisting of N — 1 binary digits, the remaining binary digit 
being interpreted as sign digit. Because it is necessary to store the sign associated 
to each partition appearing in the recursion formula, we shall represent the 
partition itself entirely in the fractional part of a word. 

First method. The partition (A) = (A1, A2, -*-,A»,) is represented by the 
number 


L = 2-1 4 Qs 4... 4p QAM rp, 


i.e., starting from the binary point, (A: — 1) zeros followed by a 1, (Az — 1) 
zeros followed by a 1, etc. The last 1, which means 2-*, is followed by zeros 
(if » < N — 1). The number ~ need not be stored separately. For example, the 
partiticn (5, 2, 2, 1,1, 1) is represented, the sign digit omitted, by the fraction 


L = .000010101111 (and then zeros). 


We use shifts to read in the components of (A), normalizations to read them 
out [2]. 

If the partitions of m are enumerated in decreasing lexicographical order, their 
representing numbers L form an increasing series. Moreover, the transition from 
one partition to the next is effected by the addition of 2-"*' to L, followed by a 
slight modification if necessary. 

This method is still used in our character program in the routine for passing 
from one partition to the next one, when several characters are to be calculated 
in succession. 

Second method. This time we write the partition 


(A) = Ay™ Ag+ > + Ag ™ 


where \; > Ax > +++ > Aw > O and all a; > 0, and put Ay — Air = 3; (Z < m), 
Am = bm. The positive integers a; and b; (¢ = 1,2, ---, m) define the partition 
completely. The partition is represented by a number L’ constructed in the 
following manner. Denoting by ¢ an arbitrary, non-negative integer (usually, 
c = 0), yet with the restriction c + a; + db; +--++ Gn + bn < N — 1, there are 
in L’, after the binary point, c ones (yet usually, c = 0), and then successively a; 
zeros, b; ones, a: zeros, b, ones, etc. Eventually, remaining space after the last 
bm ones is filled by zeros. 
Again, we read in and out with the aid of shifts and normalizations. 
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Each non-vanishing term in the recursion formula is due to the existence of 
a so-called “hook [5] of length h” in the Young diagram, h being the actuated 
component of the class partition. Such hooks are most easily recognized in the 
number L’, namely as a pair of a 0 and a 1, spaced h binary places from each other, 
the 0 being the most significant digit in that pair. The sign of the corresponding 
term of the formula is (—1)*, where z is the number of zeros between the 0 and 
the 1 which were said to form the pair. The removal of the hook is an operation 
to be performed before the recursion formula is used again. It is equivalent to 
the interchange of the two digits of the mentioned pair. 

These operations, namely the searching for hooks, the determination of signs, 
and the removal of hooks, form the essential contents of the calculation scheme. 
We note that the removal of a hook will never cause an increase of the number of 
digits of L’. Finally, the conjugate partition (X) is very easily found, being repre- 
sented by the ‘‘complement on 1” of the number obtained by reading L’ in the 
reverse order. BESK has an instruction for shifting the multiplier register into 
the accumulator in reversed order. 

Third method. The methods just described imply that an integer, ¢, will be 
represented by ¢ digits. The fact that each digit is capable of two values is used 
only to realize the boundaries between adjacent integers. A more efficient nota- 
tion is obtained as follows. The integer ¢, which is to be packed together with other 
integers, is first written in binary notation. The number, r, of requisite digits is 
then represented by 7 digits, all ones or all zeros. The two parts of the representa- 
tion thus obtained are most conveniently placed in the different halves of one 
and the same whole word, r being placed in the half word containing the binary 
point. Thus the normalization, being guided by r, gives us the value of ¢ in binary 
notation, unpacked. 

It should be observed that such symbols as 1, 01, 001, etc., can be used here to 
represent different values of ¢. In this way, the integers ¢ = 1, 2, 3, 4, 5, 6, 7, etc., 
will be represented by the symbols 0, 1, 00, 01, 10, 11, 000, etc., respectively, 
the corresponding values of r being 1, 1, 2, 2, 2, 2, 3, etc. Therefore, if r is de- 
termined so that 


a2—-1s%t< 2 —1, 


the number ¢ will be written in 7 binary digits, constituting a number usually 
read ast +1 — 2". 

This “half-binary” notation could be used for packing the integers c, a, 
bi, «++, Gm, bm, appearing in the second method. If the numbers of digits in the 
binary notations of these integers are , a1, 81, ---, @m, Bm, respectively, we write 
y ones, a; zeros, 6; ones, efc., in the binary point half word, and the binary nota- 
tions of c, @,, b:, etc., in the remaining half word. This gives us a notation for 
partitions requiring a number of digits roughly proportional to vn, as formerly 
indicated. 

Of course, the method could be iterated (writing r binarily, etc.), but this would 
not give a very efficient notation for partitions. 

Other applications. As an example, we indicate how the methods can be 
applied to the storing of the factorization of an integer, m, into powers of primes. 
Let p; (¢ = 1, 2, «++, m) be the primes dividing n, p;: > p2 >---> Pn > 2, and 
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let a; be the exponent denoting the highest power of #; still dividing n, a; > 0. 
Further, let the function f(x) mean the number of all existing primes <x. Putting 
S (bs) — f (dias) = 0: (4 < m), f(Pm) = bm, We can use the second or third method 


above to store the numbers a, bi, --+, @m, bm. These determine the factorization 
compietely. 

Stic Com&ét 
Matematikmaskinnamnden 


(The Swedish Board for Computing Machinery) 
Box 6131, Stockholm 6, Sweden 
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Determination of Steiner Triple Systems 
of Order 15 


1. Introduction. A Steiner Triple System is a set of subsets of m marks such 
that: a) Each subset contains three distinct marks, b) each pair of distinct 
marks appears in one and only one subset. It is easy to see that m = 1 or 3 (mod 6) 
[1] and that the number of subsets is m(m — 1)/6. For n = 3, 7, 9 there is one 
and essentially only one such set, all others being obtained by permutation of the 
marks. For n = 13, there are two distinct systems. 

Using SWAC for the computations we have obtained the non-isomorphic 
systems of order 15. During the course of the calculations, Prof. L. J. Paice 
called our attention to a series of papers [2] which culminate in the listing by 
F. N. Coe of a set of 80 systems of order 15 with proof that these are distinct 
and an argument tending to show that the listing is exhaustive. Since following 
the argument would entail tremendous hand calculation and since the number 
80 was below our expectation so that it seemed quite possible that there had been 
errors or omissions in the determination, we decided to continue. It turns out 
that these 80 are indeed complete and that Cole’s computations were either cor- 
rect or at any rate did not contain any errors which led to the omission of a 
valid system. In any case these results testify to Cole’s reputation for brilliant 
computation. 

2. The computation—General considerations. There are two stages in the 
calculation, the formation of the systems and the reduction of isomorphisms. 
That these must be carried out simultaneously is seen from the fact that if all 
systems on 15 marks which arise from the 80 distinct systems are counted, the 
total exceeds 6.10" (actually (14521/315)15!). 
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In practice the reduction was carried out at three places. First the total num- 
ber of systems to be computed by SWAC was reduced to the order of 10* by giving 
the machine partial systems to complete in such a way that all distinct systems 
were certain to be included in the final result. Second, in the later stages of the 
computation it was possible to inhibit output of computed systems by a simple 
test for isomorphism to previously determined systems which could be added to 
the completion routine. Finally the output systems had to be tested for permuta- 
tions taking them into other outputs. 

3. Initial data. It is clear that, given any Steiner Triple System on fifteen 
marks, by renaming we may assume that the triples (1, 2, 3)(1, 4, 5)(1, 6, 7) 
(1, 8, 9)(1, 10, 11) (1, 12, 13) (1, 14, 15) (2, 4, 6) are present. It is then easy to see 
that there are essentially only four ways to complete the triples containing 2. 
These are: 


A: (2,5, 7)(2, 8, 10)(2, 9, 11)(2, 12, 14)(2, 13, 15) 
B: (2,5, 7)(2, 8, 15)(2, 9, 10) (2, 11, 12) (2, 13, 14) 
C: (2,5, 8)(2, 7, 9)(2, 10, 12) (2, 11, 14) (2, 13, 15) 
D: (2,5, 8)(2, 7, 10)(2, 9, 12)(2, 11, 14) (2, 13, 15). 


Cole calls these respective arrangements, triple tetrad, single tetrad, hexad, 
and duodecad, since the first has three sets of four triples in the form: (1, a, 8) 
(1, y, 6) (2, a, y)(2, B, 6). The second has only one such tetrad and the third has 
two sets of the form (1, a, 8)(1, vy, 6) (1, e, ¢)(2, a, y) (2, 8, €)(2, 6, ). D has no 
combination of interlocking triples short of the full subsystem with the obvious 
omission of (1, 2, 3). 

Again, it is apparent that, in case A, we may assume the presence of either 
(3,4, 7) or (3, 4, 8). For suppose the triple (3, 4, k), k > 8; there is a set of four 
triples (1, &, 7) (1, m, 2) (2, k, m) (2, m, j) ; the permutation (8, k) (9, 7) (10, m) (11, ) 
will produce the triple (348) leaving the remainder invariant. Further, in the case 
(347) we may assume the remaining triples including 3 to be (3, 5, 6), (3, 8, 11) 
(3, 9, 10) (3, 12, 15) (3, 13, 14) for if, say, (3,5,8) occurred, the permutation 
(45) (67) will leave the triples on 1 and 2 invariant while introducing (348). 

In a similar fashion the cases considered under B, C, and D, may be reduced. 
In this reduction the automorphism group on the partial system plays an im- 
portant role. Finally, we may consider the sets in a hierarchy. Under B, no 
system which includes triple tetrads on the triples using any pair of marks (a, 5) 
to replace (1, 2) in the definition of triple tetrad above need be considered since 
a permutation of the form (1, a)(2, b)--- would reduce it to one of the A systems. 
Under C no system with a tetrad comes into consideration. Under D we can rule 
out any system containing either a hexad or a tetrad. This consideration actually 
serves to exclude D completely. 

The length to which this procedure is carried must be balanced between several 
conflicting requirements: 


1) The machine should produce as few systems as possible. 
2) There should be as few starts as possible since each new subsystem demands 
neW feeding of initial data which is time-consuming to prepare. 
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3) The initial hand computation must not be carried to a point where there 
is danger of losing completeness by errors arising from its complexity. 


The requirement 3) is less important than it appears at first glance since it is far 
easier to neglect existing automorphisms (which would create over-production) 
than to introduce non-existing ones. At any rate 2) and 3) together determined 
that no production of initial data beyond triples on 3 was considered. The initial 
starts, then, included from 14 to 19 fixed triples. For example, in the case A: 
(348) the following additional triples were used: 1. (3, 5, 6), 2. (3, 5, 9)(3, 12, 15) 
(3, 13, 14), 3. (3, 5, 10) (3, 12, 15) (3, 13, 14), 4. (3, 6, 12). On start 1, the machine 
was stopped after (3, 7, 9), (3, 7, 10), (3, 7, 11), had been considered. 

4. Exhaustive production of Steiner Triple Systems. As remarked the compu- 
tations were performed on SWAC. This is a binary machine using 36 digits plus 
sign in a word. The high speed memory consists of 256 words and the medium 
speed memory was a drum holding 4096 words. These sizes limited the choice of 
routines and the methods of storage. 

The machine is required to form all systems using the given first triples. 
Rather than a list of triples the machine actually stores and works with an inci- 
dence or usage table. Specifically, if we designate the first four 8 digit portions 
of a SWAC word by a, 8, y, 6, each mark is stored in a of a corresponding word; 
B holds the number of marks with which the given mark has appeared at that stage 
of the computation; the ith position of the first fifteen digits of y and 6 holds 
a 1 if the given mark has been used with 7 and a 0 otherwise. Thus the cell corre- 
sponding to the mark 8 in the entry A: (348) case, initially read in decimal: 
a = 008, 8 = 007, y = 241, 6 = 192. Note that an incidence of a mark with 
itself is always counted so as to instruct the machine -unambiguously that a 
mark i is available for use in a triple with a mark j if and only if a 0 appears in 
the ith place of (7, 6) in the word corresponding to 7. With this and the binary 
representations of decimal 241 (1111, 0001) and decimal 192 (1100, 0000), it is 
immediately seen that the triples (1, 8, 9)(2, 8, 10) (3, 4, 8) generate the line. 

The basic code may be outlined as follows: 


Possible Successor 
Pre- 

Routine decessors Operations Performed Criterion Yes No 
A F Find least a whose8 <15 Are all 6 > 15? B D 
B A Output (or go to supple- No choice except C C 

mentary routine for possibility of in- 
search for earlier iso- serting supple- 
morphs) mentary routine 
C  Borsupple- Return to next earlier No choice D D 
mentary value of a 
routine 
D AorC Determine marks avail- No choice E E 


able for use with pres- 
ent choice of a value 
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Possible Successor 
Pre- 
Routine decessors Operations Performed Criterion Yes No 
E D Permute to form new Isanewsubsystem F G 
subsystem on mark de- found ? 
noted by present value 
of a 
F E Store new subsystem with No choice A A 
other completed ones 
G E Haveall permutationsfor Have all permuta~- H E 
this a been tried? tions for this a 
been tried? 
H G Recall next earlier a and _ Is an earlier a ae 
its list available? 
I H Halt (end of calculation) 
A 
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The permutation scheme of routine E is a modification of a Paige-Tompkins 
routine [3] adapted to the fact that order within a triple is of no importance and 
thus there are only (m — 1)(m — 3)---1 permutations on m marks to be tried. 
Both the permutation scheme and the complete code can be illustrated by the 
example A: (348) (356). Here 1 and 2 are complete, their 6 is 15. For 3, 4, 5, 6, 8, 
B = 7. For 7, 9, 10, 11, 12, 13, 14, 15, 8 = 5. Thus 3 is selected as the least number 
still incomplete with largest possible 8. The numbers 7, 9; 10, 11; 12, 13; 14, 15; 
are available for use with 3, the semicolons indicating the initial breakdown into 
pairs which could form triples with 3. Permutation occurs as far to the left as 
possible, either because of impossible pairings or fully utilized sub-systems. Here 
(10, 11) is the first impossible one. The next considered is 7, 9; 10, 12; 13, 14; 
15, 11. This is acceptable and the lines of the table are stored in this order on a 
drum channel; the permutation tallies are also stored for future reference. The 
incidence tables are then adjusted in the high speed memory to account for this 
sub-system and we proceed. This time 4 is selected with its list of 7, 9; 10, 11; 
12, 13; 14, 15. The successive trials from this point to form a new sub-system 
are: 1) 7, 10; 11, 12; 13, 14; 15,9; 2) 7, 10; 11, 12; 13, 15; 9, 14; 3) 7, 10; 11, 12; 
13, 9; 14, 15; 4) 7, 10; 11, 12; 14, 15; 9, 12; 5) 7, 10; 11, 13; 14, 9; 12, 15. The 
italics indicate the (first) impossible pairing. Of course, after step 3) the order 
returns momentarily to 1) while the tallies are adjusted to indicate completion 
of this level and to start the permutation at the next level in the prescribed pattern. 

After storage, 7 is called with 8, 11; 12, 13; 14, 15 associated. This time no 
acceptable sub-system can be found so 4 is recalled and begins permuting where 
it left off. It happens that this entry on 3 cannot be completed ; none of the twenty 
apparently acceptable sub-systems on 4 associated with it lead to a complete 
system. Not until 3; 7, 11; 12, 15; 9, 10; 13, 14 does a completable entry on 3 arise. 

To output, the a entries from the successive storage drum channels are col- 
lected in a compact arrangement and printed or punched. The code described 
here was written by LAVERNE E. RICKARD with suggestions from C. B. TompKIns 
and the authors. 

5. Reduction to distinct systems. For the A and B levels the problem was run 
with no additional routine. Using a tabulation of the output, those systems which 
were obviously isomorphic to earlier systems but had been generated because it 
was simpler to let the machine run than to cut them out by increasing the number 
of starts (requirement 2 above) were discarded by hand. 

In the C computation, a routine was added before output. This routine looked 
at each completed system checking possible tetrad linkages. If one was found 
the output was by-passed and the next system was computed. This resulted in 
the production of only two C-systems which were found to be isomorphic on 
hand-checking. 

The problem of removing the non-obvious isomorphisms still remained. At 
this point some 5000 systems were left in the running. It was manifestly im- 
possible to try any considerable proportion of the 15! permutations on the possible 
pairs of isomorphic systems. Accordingly, using a modification of a procedure of 
Louise D. CuMMINGs contained in the cited memoir, a set of invariants was tabu- 
lated for each system. 

Given a triple system, the sub-system consisting of the triples on any mark 
possesses one of the linkage types A, B, C, D, with the sub-system on each of 
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the 14 other marks. To each mark, then, can be correlated a vector (a, b, c, d) 
where a + 6 +c +d = 14; a is the number of A linkages with the given mark, 
b the number of B linkages, etc. Thus associated with each system is a 15 X 4 
matrix and if two systems are isomorphic their matrices are row permutations 
of each other. Further, if system S; is carried into system S, by the permutation 
be Hl rity +) , the vector corresponding to j of S; must be equal to the vector 
corresponding to i; of Sp. A routine was written for SWAC which computed these 
vectors (and their sum for sorting purposes) for each system. When a system 
was input, a table was first formed which stored in the (7, 7) position the number 
k appearing in the triple (i, j,k) of the system. Zero was stored in the (#, 4) 
positions. Since 4-15 = 60 digits were required for any line of the table, double 
precision computation was necessary. The main part of the routine then took 
each ordered pair (i, 7) and found the coincidences. If there were three tetrads 
a tally was added in the a section of the ith output word, etc. In working with the 
B systems a rejection step was added whenever a triple tetrad occurred. 

There now remained 1912 A-systems and 538 B-systems which were sorted 
into stacks according to the equivalent matrices. The final routine used a repre- 
sentative of each of the stacks and attempted to find a permutation taking the 
successive systems into this archetype; i.e., choices were made for 7, and i; from 
the possible numbers with correct associated vectors. This determined 7; from 
(1, 2, 3) which was checked for vector equality. If it did not check, a new choice 
was made for 72; if it did check, a choice for 14 was made which determined 7; 
from (1, 4, 5), i¢ from (2, 4, 6), % from (34k), 7, from (35m), etc. Each of these 
was checked for vector equality and for consistency. (E.g., if m = 6, ig = im.) 
The permutation found was then output. It turned out that the invariants were 
complete in the sense that all systems with equivalent matrices were isomorphic. 

6. Time factors and controls. Systems were produced at rates varying from 
20 per minute down to three per minute depending on the relative “richness” of 
the initial starts. The linkage computation required about 2 seconds per system. 
The final permutation routine required an average of 4 seconds per system in- 
cluding typewriter output of the permutation. 

All routines used were essentially self-checking. In the exhaustive production 
routine the possibility of a consistent error in all the usage tables is negligible. 
An error will show up as “‘i is available for use with 7 but 7 is not available with 7.” 
When this happened the machine went into a short loop and could only be re- 
started by feeding consistent data. Several hand computational checks were 
made to insure that it was indeed producing the “‘next system.” It may be noted 
that every possible distinct system occurred in at least two isomorphs. Some of 
them occurred 200 times. 

The checks on the vector computation were that the matrix obtained was 
symmetric, that every vector had component sum 14, and that the matrices were 
consistent with others of the same type. 

For the permutation routine, the various requirements of consistency were 
again such as to make the possibility of finding a false permutation negligible. 
An error could, however, cause an answer of ‘‘not isomorphic’’ to appear. This 
did, in fact, happen in several cases. In each case the system was re-fed twice; 
the resulting permutations were compared and hand checked. 
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7. Conclusions. The chief source of wasted time lay in the lack of initial 
knowledge concerning the magnitude of the problem: 


a) The original output of the Rickard code was in converted decimal, designed 
for easy tabulation and hand procedures. These outputs had to be reconverted 
into binary form readily assimilated by SWAC. 

b) Systems with A-linkage should have been rejected without output during 
the B-computation as was done for tetrad linkage in the C-computation. In fact, 
the linkage vectors should have been computed simultaneously with the system. 
However, we had not thought of these matrices until these calculations were 
complete. 


The procedures described in this paper represent, to the best of our knowledge, 
the first use of large scale digital computing equipment to isolate isomorphic 
algebraic systems. In this connection, the crucial step is the devising of a reason- 
ably calculable system of invariants to reduce the gross search for mappings. By 
its nature, a computer is better adapted to calculating such invariants than to 
exhaustive search for transformations. 
The question of computing the systems of order 19 arises. Here there are seven 
possible starts on the systems using 2; these may be represented by the possible 
partitions of 16 into even numbers greater than 2 in the same way that A, B, C, D, 
in the fifteen case correspond to 4 + 4 + 4, 4 + 8, 6 + 6, 12 respectively. For 
each of these starts there are 2-10° a priori distinct permutations of (4, 5, ---, 19) 
for use with 3 compared to 10‘ on (4, 5, ---, 15). Machine procedures would be 
required to isolate the starts and the total number requiring trial might be as 
great as 10‘. The problem thus seems to be beyond the reach of present equip- 
ment. Some short runs of systems have been made on SWAC and these bear out 
this impression. Unfortunately, no reliable estimate of the total number of systems 
to be expected exists. 
The authors wish to acknowledge the assistance given them in coding and 
computation by C. B. Tompkins, L. J. Paige, Ruth B. Horgan, LaVerne E. 
Rickard, and John Selfridge. The codes used are on file at Numerical Analysis 
Research, University of California, Los Angeles. 
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A Set of Test Matrices 


1. Introduction. For the testing of the adequacy and efficiency of any com- 
puting method it is clearly very desirable to have readily available a set of reliable 
test data. The material discussed below has been found by the author to be of 
great utility for this purpose. It is surmised that it may also be of interest to 
many others engaged in the testing of new techniques for the inversion of matrices 
and the determination of characteristic roots and vectors. Furthermore, the 
results may be of interest also for purely theoretical reasons. 

The exact inverses of a sequence of ten nonsymmetric and again of ten sym- 
metric matrices of extremely poor condition, their determinants, and their largest 
and smallest characteristic roots and associated vectors as well have all been 
calculated. These are available in the UMT file. These matrices were obtained by 
a simple modification of the well known Hilbert matrix, and while they share 
many features characteristic of that matrix, they differ from it in many other 
important aspects. 

The data calculated also permit an examination of the closeness of certain 
bounds proposed by OsTrROwSKI, PARKER, eét al. 

2. The sequence of matrices. The set of matrices mentioned above is obtained 
from the sequence of “‘positive’’ matrices. 


rc 


1 1 eA 1 7} 
241 331 (n + 1)7 
3-1 4-1 (mn + 2) 
A, = (a4) = ° ° ° , 2 = 1,2,3, --:. 
n (n+ 41)> (2m — 1)*) 








Thus a;; = 1 forall 7 = 1, 2,3, ---, anda, = (¢ + 7 — 1)“ fort = 2,3, ---; 
j =1,2,3,-+-. 
The matrix A, has the following properties: 
a. The determinant d, of A, is of the form (—1)*~é,—, where 6, > 0 is an 
integer. 
b. The elements aj of A,— are integers. Further, 


~ 1 for j=1 
(s) _ J 
2 aj 0 for j #1. 


c. The characteristic root of largest absolute value is positive, that of smallest 
absolute value negative. 


The nonsymmetric matrices A, are thus finite segments of a matrix A which 
is closely related to the well known symmetric Hilbert matrix 


(2.1) B = (64), bs = @ + 7-—1)7 


for i, 7 = 1, 2,3, ---. The matrix B has been discussed in many places, e.g., by 
Todd [1], and many of its properties are now known. 
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In giving brief sketches in substantiation of some of the properties of A, we 
shall prove first that 6, satisfies the following recursion formula: 


(2.2) dos = (27, )(2") cn + 108 


for m = 1, 2,3, ---, with 6; = 1. 
Multiplying each column of d,,; by the denominator in the last row, and then 
dividing outside by the product of these factors, yields 














nm+1n+2 --- 2n+1 
nm+in+2 <--++ 2n+i 
1 2 3 n+2 
day = . . 
(m + 1)(m + 2)---(2n + 1) 
1 1 yal 1 








The following further steps, carried out in the order indicated, lead to (2.2): 
1. Subtract the last column from each of the others, making use of the 
fact that 





nt+j _mt+1_ (m—-it+1)(m—f+1) 
e+7—-1t sts (@+ j—1)\(1 +2) 


2. Remove from the first » columns the respective factors m,n — 1, ---,1 
to the outside, and reduce the determinant to the order n. 
3. Finally remove from all the rows except the first the respective factors 





a—-in*—2 . 3 
n+2'’n+3’ °2n’ 
getting 
(2.3) den = n'(n — 1)! 
° a1 = 





~ (n+ 1)(m + 2) + (Qn +i)(n +2 n+ 3)---2n™ 





--| (n — 1)! I. n P 
i (n +2) +--+ (2n)1 (2n+1)(m+1) * 





2n —2 n 
- whys :) ‘Ont Dm +1) ™ 


which proves (2.2). 
The solution of the difference formula (2.3) leads to 


[1!2!--- (m — 2)!(m — 1)!}* 





= aes n—l 
(2.4) d,, ( 1) (m — 1)!(m + 1)!(m + 2)!--- (2n — 1)! 
for m = 1, 2,3, ---. 
Consequently, 


det A, = (—1)"7 [1!2! +++ (m —1)!}! 





1!2!--- (2n — 1)! 








so 


fo 
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In [1] it was shown that for the finite segment B, of order n of Hilbert’s matrix, 


_ [alat--- (m — 11) 
det B.= “Tol (Qe — 1!’ 





so that 
det A, = (—1)"""'n det B,. 


The value of d, may be estimated quite easily using (2.2) and Stirling’s 
formula n! ~ (n/e)"(2xn)"/?; the result is 


|dn| ~~ 3°20 9g (n), 


with ¢(m) of the form Con’, Cy denoting a constant. 


TABLE 1a. Values of 5, 


= 


1 

6 

720 

15 12000 

5 33433 60000 

31 05223 67232 00000 
29541 55782 75110 09280- 104 
45669 60589 07168 10735-10" 
11430 90871 53174 41088-10” 
46206 89394 79146 91316-10” 


COCO ONAUEP WN 


— 


TABLE 1b. Measures of condition 


n d(A,) M(A,) P(A,) 

1 1 1 

2 6— 12 12-587 

3 720-4 540 354-51 

4 (15120-10?)— 17280 13090 

5 (53343-10°)— 67200-10' 45057 -10! 

6 (31052-10")— 23814-10* 15259-10* 

7 (29542 -10*)— 80681 - 104 51270-10* 

8 (45670-10?7)— 28333 -10° 17164-10° 

9 (11431 - 1087) 95447 -107 57364-107 
10 (46207 - 1047) 33640 - 10° 19158-10° 

A listing of the values 6, for m = 1, 2, ---, 10, rounded to 20 significant figures, 


is given in Table 1a. The extremely rapid increase of 6, is worthy of note; it 
serves to underline the poor ‘‘condition” of the A,, as measured in terms of d,. 
As will be shown later other, more suitable measures of condition, 


M(A) = 1n:-||A|| |A“,  ||A|] = max |ay| 
i,3 
P(A) = |Amax(A)/Amin(A)|, 
are indicative of the same phenomenon. 
Next a few remarks about the inverses A,~'. All the elements of A,~' are 








: 
\ 
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integers. In fact, it can be shown that for i = 1, 2, ---, we have 
() — ¢_y\ni(*@tt—1\(n 
of = (-4(*T471)(4) 
while for 1 = 2,3, --- and j = 1, 2, ---, ”, we have 
(nm) = aS ib, . n+i-—1 n+j 
oft = (Dose (™ FSG ETH) 
where 


Rijn. = we ert ‘Vi 


Let us proceed to the characteristic roots \; and characteristic vectors x; of An, 
arranging the \; according to increasing absolute value: |\i| <|A2| <---<|Aal; 
we shall set A1 = Am, An = Am. Let the vector x; have the components x;;, 
j =1,2,---,n. 

Since 


det [A _ AT] = (—1)"A[A"" = (2; a:;)A" +.:- -] + (—1)"-6,-1, 


it follows immediately that there is always at least one positive characteristic 
root. It seems that there is only one positive characteristic root, the others being 
negative. 

The existence of a positive characteristic root also follows from the result of 
FROBENIUS [2] that the root \y» of largest absolute value of any positive non- 
decomposable matrix is positive and simple, and may be associated with a proper 
vector having positive components only, the only vector possible that possesses 
this property. ‘ 

Estimates for the \; are readily calculated from the following theorem due to 
PARKER [3]: 


Let U = (uj;) be an arbitrary matrix of order n. Put 


g = 0S [wil P=> | wis|, = | wss| 5 
= ne hii 
S; = Pi+ lui —é|, S = max;S;; T; = Q; + |uj —&|; T = max; Tj. 
Then for every characteristic root \ of U 
(2.5) |X — | < min (S, T). 


Furthermore, the number & may be replaced by any other number. 


Applying (2.5) to the matrix A,, for = 0, it is seen that the maximum column 
sum T is a bound: 


(2.6) Au™ < log (n +1) +7, 


with y = 0.577--- denoting Euler’s constant. 
Bounds for the components of xy» may be obtained by means of Ostrowski’s 
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theorem [4]: 
Let A, = (aij) be a matrix of order n of positive elements. Put 


R; = > ay, + = 1,2, ---, 2; R = max R;, _ r = min R;: 


j=i i é 


ky 


ll 


min aii, ke = min aj;; T= C(r — ky)/(R - kx) }. 

i x7 
Further, let the positive components xm, of the characteristic vector xy belonging to 
Am be normalized so that 


XmMn S Xmjn-1 S++ S XmM2 S Xm = 1. 
Then 
ko/(R — r + ke) < XmMa S o1 


provided R # r. 
In our case 


R=n, ky = (Qn—1)~, ke = (Qn —2)7, r= (n+ j-1)-. 


= 


It is found that 
(2.7) (2n?)"[1 + O(n) ] < xmun < [log 2/n}*[1 + O(n) ]. 


The characteristic roots of least absolute value X,, are frequently also of in- 
terest. If the A; are distinct, then 


n 


A,2 = pid (1/Ag™) xe «7 
i=1 

where the x,“ denote the characteristic column vectors of A, belonging to A;, 
and the r;™ are the characteristic row vectors of the transposed matrix A,” 
belonging to A,. The smallest root A,“ may then be found as the dominant 
root uw™ of A,-". 

If it is further known that the roots A of the matrix U = (u,;) are real, then 
(2.5) permits the conclusion 


A > — — min (S, T); 
if, in addition, the particular root \ = uy is negative, then clearly 
(2.8) uu < [é — min (S, T) -. 
Another upper bound for A,, is provided by the fact that for definite matrices 


(n||A ||) < JA“ < |An|—. 
Thus 


(2.9) [Am] < [|A-". 


3. The inverses and determinants. The first ten inversions were carried out by 
the method of partitioning [5]. A tabulation of the inverse Ag computed by 
this method is given in Table 2. 
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TABLE 2. The inverse of Ag 


—- 6 630 — 6720 22680 — 30240 13860 
105 — 7350 88200 — 3 17520 4 41000 — 2 07900 
Agta | 560 29400 —3 76320 14 11200 —20 16000 9 70200 


— 1260 44100 —6 04800 23 81400 —35 28000 17 46360 


1260 —52920 7 05600 —27 21600 39 69000 —19 40400 
462 —13860 1 94040 — 7 76160 11 64240 — 5 45360 


Certain interesting properties of these inverses have already been pointed out 
before. The following features of the A,~ are also noteworthy: 


a. It is seen that 


( 


(—1)*— for i > 1, 7 for m even 
sgn af = j 


21, 
(—1)*— for z > 1, 7 > 2, for n odd 
(—1)** fori > 1, 7 = 1, for m odd. 
b. Aa] = offs. 241 or af"h4, with k +1 <n for n > 3. 


Now 
An] > la | ~ C2. 


Inspection of the ||A,~|| suggests that 
(3.1) An || ~ C.2% with C, ~ 4-10-%. 

With these estimates it is seen that the measures M(A,) are approximately 
(3.2) M(A,) ~ Cn2™, 


so that the order of ill-conditioning is about that of Hilbert’s matrix. 

The exact values of the M(A,) are shown in Table 1b. 

It has been stated [6] that “average’’ matrices have M-condition numbers of 
the order n/? log n, and P-condition numbers of the order n. In the light of such 
orders for condition numbers the matrices A, must indeed be called extremely 
pathological. 

There is a large number of techniques in matrix algebra which are specifically 
restricted to symmetric positive matrices H; for the use of such procedures there 
have been calculated the matrices 


Ay = An™An. 
From the reciprocals of A, those of H, are easily calculated as 
(3.3) HH, = A,(A,)?. 


The inverses H,,— thus also consist of integer elements only. It is known (TAUSSKY 
[8]) that the condition (in several senses) of AA’ is worse than that of A. 
Obviously in Ais: = (h§*”), 


(3.4) bf? =14+(64+)G4+)P4+164+29)G4+29PRh+-: 


+(E+h)G+4)P 
=hP +(6+H)G+O)P 
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for all i, 7 < k, k = 1,2, ---, with hf? = 1. Consequently, 


Ww =14vH- 5 G+, 
jak+i 
where W(x) = I'’(x + 1)/I'(x + 1). Further, for fixed k, 
(3.5) max hf? =A? =14+W1)-E51+)%=8/6-5 G4+H-. 
4.7 =k p=k 


A listing of H is provided in Table 3. 


TABLE 3. The matrix Hz 





( 5369 19 431 1867. 12107 1 62997) 
3600 14 336 1512 10080 1 38600 
2 22581 29 169 1241 17911 
1 76400 24 144 1080 15840 
8 22949 41 403 6563 
7 05600 36 360 5940 
Hy = 
70 88041 11 359 
63 50400 10 330 
68 97529 71 
63 50400 66 
8196 07009 
L 7683 98400 | 








4. Characteristic roots and vectors. Since there is a unique characteristic root 
Amu™ of largest absolute value the power method [7] may be used to isolate 
Aw™, while simultaneously obtaining the associated characteristic vector xy™ : 


Agxe?, -” Af? xf”, k= 1, 2, 3, _" 


lim k™ = Av”, lima” = xu, 
kon kw 


(4.1) 


TABLE 4. Characteristic roots 


n Am™ An” 

1 1.000 000 

2 1.448 403 —.115 0693 

3 1.707 105 —.481 5399-10 
4 1.886 632 —.144 1324-10 
> 2.022 999 — 448 9833-10-5 
6 2.132 376 —.139 7499-10-* 
7 2.223 362 — .433 6577-10-8 
8 2.301 055 —.134 0623-10 
9 2.368 717 —.412 9309-10-4 
10 2.428 554 —.126 7649-10-" 
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provided the initial arbitrary vector xo is not orthogonal to xy‘. The rate of 
convergence increases with the ratio between Ay and the next characteristic 
root of A,. 

The absolutely largest elements of A, occur in the first row; it is therefore 
advisable to normalize the xy” by keeping its first element equal to unity. In 
carrying out the computation a considerable saving in labor may be further 
achieved by utilizing the values of xy” obtained at the mth step to start the next 
step with a vector xo'*t) whose first m components were identical with xy. 
Proceeding in this manner there was obtained in seven or less iterations an agree- 
ment to six decimal places in both characteristic roots and vectors. 

The characteristic roots Ay“ are shown in the middle column of Table 4; 
as indicated by (2.6) they actually do not grow faster than log n. 

The characteristic vectors xy™ belonging to the As are listed in Table 5a. 
Their smallest components x{P, seem to decrease as n—. The bounds given in 
(2.7) thus seem a bit wide. 

The characteristic roots \m‘” = ua of least absolute value were also calcu- 
lated, together with their characteristic vectors x,“, now normalized to have 
their last component equal to unity. The values of X,,“ are shown in the third 
column of Table 4. Here the inequality (2.8), with — = 0, is found to restrict the 
roots fairly well, the upper bounds—min (S™, T™)— being, at least for m < 10, 
less than 30 per cent too high. 

For purposes of comparison the resulting condition numbers P(A,) were also 
computed ; they are exhibited in the last column of Table 1b. These suggest the 
following rough estimate: 


(4.2) P(A,n) ~ C32" log n, Cs ~ 8-10-%, 
which, with (3.2), gives, 

(4.3) P(A,) ~ (2/n) log nM(A,); 
Todd [1] found for Hilbert’s matrix B 


anM(B,) < P(B,) < +M(B,). 


Mark LOTKIN 
AVCO Manufacturing Corporation 
Stratford, Connecticut 
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Some Comments on a NORC Computation of 


Among the problems suggested as demonstration routines for the NORC [1] 
was the calculation of x to a large number of digits. In 1949, x was calculated to 
2035D on the ENIAC [2]. The present computation on the NORC, which was 
carried out by the authors at the Watson Scientific Computing Laboratory, pro- 
duced 3089D. This limit was chosen since the entire calculation could be contained 
in the NORC’s high speed memory (2000 locations). The program was designed 
to produce any number of digits up to this limit. More extensive programming 
involving use of the magnetic tapes would have increased the number of digits 
by several orders of magnitude. 

The formula used was the same as that employed for the ENIAC computation : 


a/4 = 4arc tan 1/5 — arc tan 1/239 


in conjunction with the GREGORY series 
arctanx = >> (—1)"(2m + 1)—1x**", 
n=0 


The form of the expression used for the NORC computation was: 


100(0.2)2"+* — ofr | 
2n +1 , 





~ Di (| 


All intermediate factors were truncated without rounding to 3093D. For 
n = 0, the value of 100(0.2)?** was set at 0.8; the value of this factor at each 
succeeding m was computed by multiplying the value at m — 1 by .04, with 
truncation of the product to 3093D. For nm = 0, the value of (1/239)?**+" was 
computed with truncation of the quotient to 3093D; the value of this factor at 
each succeeding was computed by dividing the value at m — 1 by 2392, with 
truncation of the quotient to 3093D. For each n, the tens complement of this 
second factor was taken and added to the first factor; the sum was divided by 
2n + 1, with truncation to 3093D. This result or its tens complement, according 


n—1 


to the value of m, was added to >>. This process was continued until the first 
0 


3093D of the second factor were zeros; computations involving this factor were 
then bypassed. The process terminated when the first 3093D of the first factor 
were zeros. A count was kept of the number of terms required for the calculation; 
at termination, m equaled 2214. 

A second run, with the same program, of the calculation was made about 
two months after the first run. A comparison of the two results showed agreement. 
The first 2035D of the ENIAC and NORC results agree. Since the programming 
for the computation on the NORC is completely general for any number of digits, 
the correctness of the digits up to 3089D is strongly indicated. Further assurance 
of the correctness of the result is afforded by the NORC’s inherent checking 
system [3]. Note that the only expected error results from truncation of the 
intermediate calculations; this error is a function of m and could affect the last 
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3.14159 
58209 
82148 


48111 7 


44288 
45648 
72458 
78925 
33057 
07446 
98336 
60943 
00056 
14684 
42019 
51870 
50244 
71010 
59825 
18577 
38095 
03530 
55748 
81754 
85836 
94482 
93313 
25338 
67823 
55706 
32116 
63698 
81647 
16136 
45477 
56887 
82796 
73929 
06744 
46776 
94657 
26719 
49625 
57098 
68683 
84598 
43904 
41694 
01684 
64565 
15076 
22874 
90097 
22658 
54285 
21359 
03742 
86875 
81911 
98391 
56794 
66389 
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37787 
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89793 
59230 
32823 
84102 
66593 
34603 
06315 
01133 
57595 
62749 
44065 
05392 
45263 
22495 
21290 
49999 
34690 
78387 
28755 
17122 
10654 
68995 
45415 
49393 
63707 
77472 
89891 
35587 
63600 
85054 
87202 
54252 
61452 
52552 
86251 
04946 
81454 
08412 
20391 
73962 
95126 
84826 
49399 
41059 
27741 
44695 
13654 
58484 
52267 
35488 
94510 
60101 
67598 
75640 
95265 
23144 
05785 
06430 
95206 
95618 
95146 
08303 
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23846 
78164 
06647 
70193 
34461 
48610 
58817 
05305 
91953 
56735 
66430 
17176 
56082 
34301 
21960 
99837 
83026 
52886 
46873 
68066 
85863 
77362 
06959 
19255 
66010 
68471 
52104 
64024 
93417 
94588 
75596 
78625 
49192 
13347 
89835 
01653 
10095 
84886 
94945 
41389 
94683 
01476 
65143 
78859 
55991 
84865 
97627 
06353 
46767 
62305 
96596 
50330 
52613 
14270 
86782 
29524 
39062 
21845 
14196 
14675 
55022 
90697 


26433 
06286 
09384 
85211 
28475 
45432 
48815 
48820 
09218 
18857 
86021 
29317 
77857 
46549 
86403 
29780 
42522 
58753 
11595 
13001 
27886 
25994 
50829 
06040 
47101 
04047 
75216 
74964 
21641 
58692 
02364 
51818 
17321 
57418 
69485 
46680 
38837 
26945 
04712 
08658 
98352 
99090 
14298 
59772 
85592 
38367 
80797 
42207 
88952 
77456 
09402 
86179 
65549 
47755 
10511 
84937 
19838 
31910 
63428 
14269 
52316 
92077 


83279 
20899 
46095 
05559 
64823 
66482 
20920 
46652 
61173 
52724 
39494 
67523 
71342 
58537 
44181 
49951 
30825 
32083 
62863 
92787 
59361 
13891 
53311 
09277 
81942 
53464 
20569 
73263 
21992 
69956 
80665 
41757 
72147 
49468 
56209 
49886 
86360 
60424 
37137 
32645 
59570 
26401 
09190 
97549 
52459 
36222 
71569 
22258 
52138 
49803 
52288 
28680 
78189 
51323 
41354 
18711 
74478 
48481 
75444 
12397 
03881 
34672 


50288 
86280 
50582 
64462 
37867 
13393 
96282 
13841 
81932 
89122 
63952 
84674 
75778 
10507 
59813 
05973 
33446 
81420 
88235 
66111 
53381 
24972 
68617 
01671 
95559 
62080 
66024 
91419 
45863 
90927 
49911 
46728 
72350 
43852 
92192 
27232 
95068 
19652 
86960 
99581 
98258 
36394 
65925 
89301 
53959 
62609 
14359 
28488 
52254 
55936 
79710 
92087 
31297 
79641 
73573 
01457 
08478 
00537 
06437 
48940 
93014 
21825 


41971 69399 
34825 34211 
23172 53594 
29489 54930 
83165 27120 
60726 02491 
92540 91715 
46951 94151 
61179 31051 
79381 83011 
24737 19070 
81846 76694 
96091 73637 
92279 68925 
62977 47713 
17328 16096 
85035 26193 
61717 76691 
37875 93751 
95909 21642 
82796 82303 
17752 83479 
27855 88907 
13900 98488 
61989 46767 
46684 25906 
05803 81501 
92726 04269 
15030 28618 
21079 75093 
98818 34797 
90977 77279 
14144 19735 
33239 07394 
22184 27255 
79178 60857 
00642 25125 
85022 21066 
95636 43719 
33904 78027 
22620 52248 
43745 53050 
09372 21696 
61753 92846 
43104 99725 
91246 08051 
97700 12961 
64815 84560 
99546 66727 
34568 17432 
89314 56691 
47609 17824 
84821 68299 
45152 37462 
95231 13427 
65403 59027 
48968 33214 
06146 80674 
45123 71819 
90718 64942 
20937 62137 
6259 


37510 
70679 
08128 
38196 
19091 
41273 
36436 
16094 
18548 
94912 
21798 
05132 
17872 
89235 
09960 
31859 
11881 
47303 
95778 
01989 
01952 
13151 
50983 
24012 
83744 
94912 
93511 
92279 
29745 
02955 
75356 
38000 
68548 
14333 
02542 
84383 
20511 
11863 
17287 
59009 
94077 
68203 
46151 
81382 
24680 
24388 
60894 
28506 
82398 
41125 
36867 
93858 
89487 
34364 
16610 
99344 
45713 
91927 
21799 
31961 
85595 
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four digits of the final result. The result has been tabulated to 3089D; the final 
digit is unrounded. 

Running time for the 3093D was approximately thirteen minutes. The pro- 
gramming takes account of the number of zeros generated to the right of the 
decimal point in each factor, so that the number of operations required for each 
term in the series decreases. This leads to the following statement—if the time 
to compute z to m digits is ¢ units, then the time to produce km digits is roughly 
k*t units; this holds true as long as the calculation is contained in high-speed 
storage. 

The following table gives a count of each of the digits in z. 


(1) (2) (3) (4) (5) 
1-3090 1-2036 2037-3090 (4)/(3) 
0 269 184 85 46 
1 315 213 102 AT 
2 314 210 104 50 
3 276 191 85 45 
4 322 198 124 63 
5 326 211 115 54 
6 311 204 107 52 
7 297 200 97 49 
8 318 207 111 54 
9 342 218 124 57 
> 3090 2036 1054 52 


S. C. NICHOLSON 
J. JEENEL 


- 


Watson Scientific Computing Laboratory 
612 West 115th Street 
New York 27, New York 


1. The IBM-Naval Ordnance Research Calculator, now located at Naval Proving Ground, 
Dahlgren, Virginia. 

2. GEORGE W. REITWIESNER, “An ENIAC determination of x and e to more than 2000 decimal 
places,” MTAC, v. 4, 1950, p. 11-15. 

3. For a description of the NORC checking system, see W. J. Eckert & R. B. Jones, Faster, 
Faster, McGraw-Hill Book Company, New York, 1955, p. 98-104. 


Orthogonal Polynomials Arising in the Numerical 
Evaluation of Inverse Laplace Transforms 


Abstract. In finding f(#), the inverse LAPLACE transform of F(p), where (1) 
ctio 

f@® = (1/2mi) f : e?'F(p)dp, the function F(p) may be either known only 
c—jo 


numerically or too complicated for evaluating f(t) by Caucuy’s theorem. When 
F(p) behaves like a polynomial without a constant term, in the variable 1/p, 
along (c — j~,c + j@), one may find f(¢) numerically using new quadrature 
formulas (analogous to those employing the zeros of the LAGUERRE polynomials 
in the direct Laplace transform). Suitable choice of p; yields an n-point quadrature 
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formula that is exact when p2, is any arbitrary polynomial of the (2m)th degree 


in x =1/p without a constant term, namely: (2) a/2np) €?pon(1/p)dp 
c—jo 
= DY As™pon(1/p;). In (2), x; = 1/p; are the zeros of the orthogonal polynomials 
i=1 


p.(x) = II (x — x,) where (3) (1/277) 


ets 
i=1 c— 


‘ e?(1/p)pn(1/p)(1/p) ‘dp = 0,1 = 0, 
1,---,a#—1 and A;™ correspond to the CHRISTOFFEL numbers. The nor- 
malization P,(1/p) = (4n — 2)(4n — 6) --- 6p,(1/p), m > 2, produces all in- 
tegral coefficients. P,(1/p) is proven to be (—1)"e-*pd"(e/p")/dp*. The 
normalization factor is proved, in three different ways, to be given by (4) 


c+ja 
(1/2ej) [ e?(1/p)P,.(1/p) Pdp = 4(—1)*. Proofs are given for the recur- 


rence formula (5) (2%—3)P,(x) =[(4n—2)(2n—3)x+2 ]P._(x)+(2n—1)Pa_2(x), 
for n>3, and the differential equation (6) x*P,""(x)+(x—1)P,’ (x) —n’P, (x) =0. 
The quantities p;, 1/p;~ and A,” were computed, mostly to 6S — 8S, for 
4 = 1(1)n, n = 1(1)8. 

I. Introduction : Occurrence of inverse Laplace transforms. For a given func- 
tion of p, F(p), which is the direct Laplace transform of some unknown function 
f(t), for t > 0, one usually finds the f(t) from the following explicit expression : 


(1) f= 55” enroap. 


Formula (1) is known as the inverse Laplace transform of F(p). In (1) the quan- 
tity c is a real constant +0 that is greater than the real part of all the singular 
points of F(p). In practice c is usually positive, but c can be negative as long as 
for f(t) satisfying Dirichlet’s conditions in any finite positive interval the integral 


f e~*'f(#)dt is absolutely convergent (H. S. Carstaw and J. C. JAEGER [1 ]). 
0 


A note by the referee follows this paper and indicates relations between the 
present work and work published elsewhere. 

The examples treated in most textbooks on operational calculus and Laplace 
transforms contain such functions F(p) that their poles and branch points (and 
residues also) are obtainable without too much difficulty, and the inversion 
integral in (1) is evaluated by suitable deformation of the path of integration, 
and the use of Cauchy’s theorem. But there are countless other examples where 
F(p) might be too complicated to yield explicit information about the location 
and nature of its singularities without a prohibitive amount of labor. For instance, 
one will recall that in most textbook examples treating the solution of ordinary 
and partial differential equations by operational means, the original system of 
differential equations is transformed into a system whose solution F(p) is usually 
some known elementary function or a very extensively tabulated function of a 
simple differential equation (like a Bessel function), so that its analytic character 
and singularities are well known. But in actual practice one might not be fortunate 
enough to obtain such a comparatively simple F(p). Thus the transformed differ- 
ential equations might not yield a known function. Instead it might be amenable 
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only to solution in series or numerical integration, with an F(p) that is given as 
a tabulated function of p. Also, the solution F(p) might be given explicitly in 
closed form as a combination of integrals of such complicated analytic expressions 
that it might be easier to evaluate it for different numerical values of p than to 
find its poles, residues, and branch points. 

The purpose of this present article is to discuss the properties of a new set 
of orthogonal polynomials which can be the basis for convenient formulas for 
approximating f(t) in (1) for different positive values of t when one has an F(p) 
that is too complicated to show its analytic character, but which can be calculated 
for any p. 

All further discussion will now be for F(p) assumed to be exactly of the 


>. aca ; 
form >> —, i.€., a polynomial in 1 /p without a constant term. 
r=1 
To obtain a definite integral without a parameter ¢ in the exponential term, 
which is the ‘‘weight function,”’ let pf = u in (1), so that we obtain 


(1’) fo -— {"" er(*) du, 


2arjt J ico 


u\. ° aah . 
where F (*) is still a polynomial in 1/u, without a constant term. 


II. Use of orthogonal polynomials. At this point one may recall the application 
of the theory of orthogonal polynomials to quadrature formulas for definite 
integrals where the integrand is the product of a preassigned weight function and 
a polynomial P(t). There it is possible to employ the value of P(t) at m fixed 
irregularly spaced points ¢;, 7 = 1, 2, ---,, such that the resulting quadrature 
formula is exact when P(t) is any arbitrary polynomial of (2m — 1)-th degree. 


Thus for the direct Laplace transform of P(#), namely f . e~?*P(t)dt, which is 
0 


essentially f e~*Q(t)dt for polynomial Q(#), the points ¢; are taken equal to the 
0 


zeros of the Laguerre polynomials, which have been tabulated extensively (H. E. 
SALZER and R. ZucKER [2]). In the present case, even though we are not dealing 
with a polynomial in p, we can still solve the problem of finding a Gaussian-type 
quadrature formula for (1’) of approximately double the degree of accuracy of an 
ordinary quadrature formula based upon the same number of equally spaced 
points. 

Thus let p2.(1/p) be any arbitrary (2m)-th degree polynomial in the variable 
1/p, which vanishes at 1/p = 0. Consider m distinct points 1/p;,7 = 1,2, ---,m 
other than 1/p = 0 and construct the (m + 1)-point Lagrangian polynomial 
Pal esea (of the n** degree in 1/p), to pon(1/p), based upon the points 1/9,, 
4=1,2,---,mand1/p = 0. The (m + 1)th point 1/p = 0 is needed in order to 
provide for the property that p2,(1/p) vanishes at p = ©. We have for this 
polynomial approximation L‘* (1/p) the explicit expression 


(2) Lit) ( *) a“ rs Leow (4) Pon ( +), 
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where on the right hand side of (2) 


n+l1 n+1 
mw zie (1) 2 P(4-1)/ip(4-4), 
p bmi \P Pr kei \Pi Pe 

the JJ’ denoting the absence of k = i. In (2), a4: is ©, so that there is actually 
no (m + 1)-th term in (2) and L&{"(1/p) is not used. (The summation in (2) 
is written with m + 1 instead of m to avoid confusion with the (m — 1)-th degree 
coefficients L;™(1/p) which differ from the L**(1/p) by not having the 
factor p;/p.) 

Following the method in G. Szec6 [3], we consider the (2m)-th degree poly- 


nomial in 1/p, namely, pon(1/p) — L“*» (1/p) which vanishes at 1/p = 0, 1/,, 
4 = 1,2, ---,m, and thus has 


sa(3)=pii(5-2) 


i=1 


as a factor. Writing 


1 1 1 1 1 
4 (4) = z0(4) 42 (4) n— (), 
(4) pe > > p? > Tn-1 > 
it follows that 


(5) 1 ie J ( 1 ) d 1 4 owt) ( 1 ) q 
onan €? pon — = — e n ad 
21j Jc—ieo = p p 22j Jc—iw p P 


: <> ( ) ( ) 
-— P— ba\ — | rai | — } dp. 
+35 bi a ae Nes 
Thus if the second term in the right member of (5) always vanishes, (5) will be 


an m-point quadrature formula that is exact for any (2m)-th degree polynomial in 
1/p without a constant term, namely, 


© 1 [om (1) Eaen (2). 
2aj Je-ie p imt Pi 
where the “Christoffel numbers” A, are given by [6] 
(7) A hf corm (:) dp. 
27 J c—iwo p 


A sufficient condition for (6) to hold is obviously the “orthogonality” of 
(1/p)pa(1/p) with respect to any arbitrary p»_:(1/p), namely, 


1 ee 8 1 r\* . 
(8) xj bas o* 9. (2) (:) dp=0, 1=0,1,---,n—1. 


The necessity of (8) is also obvious from (6) by choosing 


Pon(1/p) >= (1/p) Pa(1/P)pn—i(1/P) 
where pn_1(1/) is any arbitrary polynomial in 1/p of the (m — 1)-th degree. 
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Hence the points 1/p;, now denoted by 1/p,), are the zeros of a certain set 
of orthogonal polynomials in the variable 1/p. 

The condition of orthogonality (8) is also mathematically equivalent, in 
terms of actual polynomials (by setting x = 1/p), to having a polynomial of the 
n* degree gn(x) which is orthogonal to any pa_,(x), with weight function e'/*/x, 
where the path of integration is a circle of radius 1/2c whose center is at (1/2c, 0). 

If the polynomial ,(1/p) is written as 


(4) + 4)" +5 (¢)" Pe (*) + 
p n—1 p n—2 p 1 p 0, 
the determination of b;,z = 0,1, ---, ” — 1, to satisfy the conditions of orthogo- 


nality (8), making use of 


1 et+jo eP 1 
a dp os << 
24j c—jeo prt 





is in the solution of this system of linear equations: 


r 











1 bn-1 Dn—2 bi , do 
tat a-s”6l[UCCU ee 
1 ba-1 ba—2 bi bo 
Ties a om 8 A 
1 bn-1 bn—2 by bo 
| (2n — i dee tes) CU ae 
For numerical work it is somewhat easier to solve (9’) in the form 
C1 + mbar + n(m — 1)dn-2 + --- + 2!b; + n!bo = 0 
1+ (m+ 1)dn1 + (a + 1)nbn2 + --- 
(9) 4 + (n+ 1) --+ 3b, + (n + 1) --+ 205 = 0 


1 + (2m — 1)b,_1 + (2m — 1)(2m — 2)bn2 + --- 
+ (2m —1)--- (n+ 1)b; + (2m — 1) ---+ nbo = 0. 





III. Explicit expression for orthogonal polynomials. It is convenient to nor- 
malize the polynomials p,(x), where x = 1/p, by multiplying p,(x), for > 2, 
by (4n — 2)(4n — 6) --- 6. This normalization produces polynomials with all 
coefficients integral (proven below) and it is not the usual normalization by 
multiplication by 


1 (te 1 i1\|/?_ 
Fa ‘ A (*)| ap] 


Denoting (4n — 2)(4n — 6) --- 6p,(1/p) by P,(1/p) for m > 2, and p,(1/p) by 
P,(1/p), one can avoid the labor of solving (9’) or (9) directly by showing that 
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P,(1/p) has the following more elegant definition : 


“0 r.()) crew (8) 


That (10) yields the leading coefficient of 1/p* in P,, (1/p), namely 


1, for #2 =1, 
(4n — 2)(4n — 6) --- 6, for n> 2, 


is obvious by induction. To prove the orthogonality property, or (8), it suffices 
to prove the vanishing of 


© BfTeilemns(e)licmn eo 
2aj Jem “> TP op pe Oi 


for m <n. This last expression is written as 


Sot" at ta —pPApmtn—l “(< -) = (< “)a 
(B) 2nj Ce a dp™ oP, 


and after integrating by parts m times, noting that the integrated parts always 
vanish, we have 


ead ape ae = (<)\¢ 
(C) a chen” Bee 


which by LEIBNITz’s rule is expressible as 


cere" ff ct+jo ™ eo )z eS d™- —r+n (<)< 
(D) 24j K ai dp’ (e- p “ ) dp» p™ dp 




















Application of Leibnitz’s rule a second time to 





Ppmtn—l 
= ~ (e-opmtn-t) 


in the above and cancellation of e?/p™, yields 


(—1)™(—1)™** (ctio ™ m+n —1 
«@ — zr E (-1)" sip 
WJ CIO ren r s=0 s 

qd™—T+" (<) 
x — } dp. 
dp=** p” p 
Now we integrate by parts (m — r + n) times each term of the above double 
summation. The integrated part will always vanish since it will have a factor of 
1/p to at least the first power. Furthermore, at some stage in the partial integra- 
tion of each term, that stage varying with the term, the integral part will also 
vanish if m < n. This last follows because the lowest power of p*~'~* is positive 


or zero, since s can equal at the most r which can equal at the most m < n — 1. 
Then in the integration by parts the positive or zero power p”~'-*, for each value 
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of s between 0 and 7, would always be eventually annulled because the initially 
occurring differential operator d™-*+*/dp™—"** is of orderm —r+n>n—1-—s 
even for the highest value of m — 1 — swhens = 0 (duetom —r+n>n-—1 
for every r between 0 and m). Thus (E) vanishes, which proves (10), and estab- 
lishes at the same time that this normalization yields all integral coefficients 
for P,(1/p). 

IV. Normalization factor. To obtain the normalization factor, which turns out 
to be given by 


e+jo 2 
(11) * a #i[p, (¢) dp = 3(—1)*, 


we repeat the preceding argument for m = n and now notice that in the final 
integral (E) the lowest power of »*-'-* will survive the integration by parts, 
because it is equal to 1/p. Retaining in the double summation in (E) only the 
single non-vanishing term s = r = m = n, we get 


(—1)?"+* ct+jo (™ ian ) m | qd" (< ) 
F) em fo Oty) SS (E) av. 


which is integrated by parts m times, the integrated part always vanishing, to give 


cq AEA 1),, fom ors Oe, 


2aj n —joo pen 








But (G) is 


2n — 1 1 (—1)"(2n — 1)(2n — 2)---n 
(—1ye( s ) min (2m)! 2 et <T 





which reduces to (—1)""/2n or $(—1)*, thus proving (11). 
From (10), the explicit formula for P,,(1/p) is seen to be [7 ] 


3 2n — 1 
on “(nen (-1) ( * )m 
sG)(Es)eon 
cm(i)(irsdunn 
com (temo 


pr 








- 


- 








n?(—1)! 
—{ 
+ ——-—- > + (—1)° 
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V. Recurrence formula. It is easy to obtain the recurrence relation for the 
polynomials P,(x) by employing a fundamental theorem about the existence of 
a recurrence formula connecting any three successive orthogonal polynomials 
(G. Szegé [8]), namely, 


(13) P(x) = (ax + bn) Pa—1(x) + CaPn—2(x). 


Thus a, is immediately seen to be 4n — 2. Equating constant terms in (13), one 
finds c, = 6, + 1, and after substitution into the equation derived from the 
coefficients of x, one obtains 


2 _ af 


b, = ’ Tes _% 
a ede PE 








so that the recurrence formula satisfied by P,(x) is seen to be [9] 


(14) (2m — 3)P,(x) = [(4n — 2)(2n — 3)x + 2)P,_1(x) + (2m — 1)Pa_2(x), 
for n> 3. 


From (14) and (8) only, without making use of (10), one can again find the 
normalization factor given in (11), through the following inductive argument: 
Multiply (14) by P,-2(x) and then operate with 


1 ct+j~ 1 


a Senet p—...@ 
2aj c—ja p p 


to obtain (making use of (8)): 


“Vise to — 3Gs = 3) ot pa (4)4p_.(4) ap 
2aj c—je p p/p P 


(2n — 1) i" 1 | ( 1 yf 
7B 4 P— i Post — dp. 
Qnj Jense “pl *\p)] 


Denoting the left member of (11) by F,, still making use of (8) to replace in the 
first of the above integrals (1/p)P,_2(1/p) by 


1 1 1 1 
<P +) ss Po (*), 
An—1 (G 4n — 6 ' p 


_ (4n — 2)(2n — 3) 
my 4n — 6 





one now obtains 


0 





F,-1 a (2n = 1)F, —2) 


or Fy; = — Fy-2. Since F; = — 3, (11) follows by induction. 
The normalization given in (11) can be seen in a third way, directly from the 
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explicit formula for P,,(1/) in (12). For, in view of (8), it suffices to consider only 


1 cti@ ep cw" - ') (*) 
H — —— iP, | — } dp, 
— Qnjdeio pp wa sat 





n 
or 

re | 2a —-r—1i1 all 
) ce) at 1) Glare )(m ws 
or 


2n—1 > ,m(m—1)---(m—r+1) — (2m—r—1)---(m+1)n 
w ( n Jat (-» r! (2n—r)(2n—r—1)---(n+1)n!' 





which, after cancellations, is written as 








. . ,(2n—1)(2n—2)---n nm n(m—1)---(n—r+1)(n—r)! 
(K) 2 (-}) (n—r)\n! 2n—r r! ? 
or 

° —— (2n — 1)(2n — 2)---n 
(L) 2 (- Vo (n — r)ir! 


and this, in turn, is expressible in the form 





(M) 3(-1)" 
x“ > (2m —0) (2n—1) (2m—2)- - - (2n—[r—1]) (2n—[r+1])---(2n—n) 
pow (r—0)(r—1)---(r>—[r—1]) (r—-[r+1])-- - (r—20) 


In (M), the $(—1)" is multiplied by the sum of the coefficients of the Lagrangian 
interpolation polynomial for the (m +1) points 0,1, ---,m, for the variable 
equal to 2m. But that sum is identically equal to 1, i.e., for any value, 2 or 
otherwise. Thus we obtain once more $(—1)* for the normalization. 

VI. Integral coefficients. It may be of interest to show that (14) alone, without 
any knowledge of (10), implies that P,(x) has integral coefficients. We prove 
this by noting that P,,:(x) will have integral coefficients if P,,(x), m <n, has 
integral coefficients and the following identical polynomial congruence holds 
form =n-+1: 


(15) 2Pm—i(x) + (2m — 1)Pno(x) = 0 (mod (2m — 3)). 


Now the existence of integral coefficients of P,,(x) and congruence (15) can be 
verified for the first few values of m. We then show that if (15) holds for some 
particular m = n, it holds for m = n + 1, provided P,,(x), m < nm — 1, has 








ir 
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integral coefficients, or, in other words, that 
2P,,(x) + (2n + 1) P,-1(x) =0 (mod (2n = 1)). 


This last congruence, by (14), is equivalent to 


4 2(2n — 1 
mag ht) + aaa Py_2(x) + (2n + 1)P,-1(x) = 0 (mod (2n — 1)), 
or to 
(2m — 1)? 2(2n — 1) ~~ 
2 wey g Py-i(x) + ~~? y P,-2(x) =0 (mod (2m — 1)), 


which in turn is expressible as 





(28 = 1) | (2m — 1)Pas(x) + 2P se) | io” Ged Gs =), 





2a.— 3 
or 
(an — 1) | (2 + (2m — 3)Py1(x) +o 1 — (2n - Pes) | ke 


(mod (2m — 1)). 


But under the assurptions that (15) holds form = n, and that P,,(x),m <n — 1, 
has integral coefficients, the last quantity in brackets is a polynomial with integral 
coefficients, which shows that the last congruence is satisfied identically in x. 
Thus (15) holds for m = n + 1 and P,,:(x) has integral coefficients. We proceed 
in this way to every n. There is a slight subtlety in the argument of this induction 
in the sense that the integral coefficients of P,,(x) up to m = nm — 1 only are 
needed to go from m = n tom = n + 1 in (15), but then use is made of the in- 
tegral coefficients of P,(x) in using (14) with » + 1 in place of n. 

VII. Differential equation. It is easy to show that P,(x) satisfies the differ- 
ential equation 


(16) x*P,"" (x) + (x — 1)P,’(x) — n*P,(x) = 0. 


Thus one merely expresses (12) in the form 





(12) P(x) = (—1y| 1 + > (—1)"n?(m? — 12)(m? — 2?) --- (nw? —r —1 =| 


r! 


and then observes that (12’) is equivalent to the automatically terminating 
“infinite series.” 


(12”) P.(x) = ¥ ax’, 
r=0 
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where 


@ = (—1)", and ra,= — (m?—r—1)a,1, for r>0. 


Working backwards from (12”), by equating coefficients of x*-, one sees that 
(12’’) must arise from (16). 


VIII. Explicit expressions for polynomials. Because these polynomials P, (x) 
are of fundamental importance, and their role in the inverse Laplace transform 
is comparable to the role of the Laguerre polynomials in the direct Laplace trans- 
form, their explicit expressions are given below for nm = 1(1)12: 

Pi(x) =x-1 

P.(x) = 6x? —4x+1 

P3(x) = 60x? — 36x? + 9x — 1 

Pi(x) = 840x* — 480x* + 120x? — 16x +1 

Ps(x) = 15120x5 — 8400x* + 2100x* — 300x? + 25x — 1 

Pe(x) = 3 32640x* — 1 81440x* + 45360x* — 6720x* + 630x? — 36x +1 

Pi(x) = 86 48640x7 — 46 56960x* + 11 64240x5 — 1 76400x* + 17640x* 
— 1176x? + 49x —1 


P3(x) = 2594 59200x* — 1383 78240x7 + 345 94560x* — 53 22240x5 
+ 5 54400x* — 40320x* + 2016x* — 64x + 1 


Ps(x) = 88216 12800x* — 46702 65600x* + 11675 66400x7 — 1816 21440x* 
+ 194 59440x5 — 14 96880x* + 83160x* — 3240x? + 81x — 1 


Pyo(x) = 33 52212 86400x'° — 17 64322 56000x*® + 4 41080 64000x* — 69189 
12000x? + 7567 56000x* — 605 40480x* + 36 03600x* — 1 58400x* 
+ 4950x? — 100x + 1 


P4;(x) = 1407 92940 28800x" — 737 48683 00800x'° + 184 37170 75200x* 
— 29 11132 22400x* + 3 23459 13600x’ — 26637 81120x* + 1664 
86320x5 — 79 27920x* + 2 83140x* — 7260x? + 121x — 1 


P32(x) = 64764 75253 24800x!? — 33790 30566 91200x" + 8447 57641 72800x" 
— 1340 88514 56000x*® + 150 84957 88800x* — 12 70312 24320x’ 
+ 82335 05280x* — 4151 34720x5 + 162 16200x* — 4 80480x* 
+ 10296x? — 144x + 1. 


IX. Zeros and Christoffel numbers. In the numerical table below there are 
given the values of the reciprocals of the zeros of P,(x) or p;, the zeros of 
P,,(x), or 1/p;, and the corresponding Christoffel numbers A, for m = 1(1)8. 
Use of these quantities in the quadrature formula (6) above can give theoretically 
exact accuracy for any polynomial in 1/p (with no constant term) up to the 16th 
degree. However, the fact that these tabulated values of p;, 1/p; and A,™ 
are correct to only about a unit in the last significant figure that is given, must 
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be taken into account in any practical example where some upper bound for the 
error should be estimated. 


HERBERT E. SALZER 
Ordnance Corps 
Diamond Ordnance Fuze Laboratories 
Washington, D. C. 
NOTE BY REFEREE 


The function F(p) is subject to certain restrictions because it is a Laplace transform. In order 
for F(p) to be the Laplace transform of the function {(¢) given by (1), it is sufficient that F(p) have 
the form (cf. G. Dortscn [5]): 

F(p) = a/p + Fi(p)/p™, 


where 6 > 0, a is a constant, and F;() is analytic and bounded in the half plane —y >. 
We assume that this condition is satisfied. Whether this condition is also sufficient for the con- 
vergence of the #-point quadrature formula to the true value of f(¢) in (1), when m tends to infinity, 
has not been determined. The author makes use here of the fact that the convergence occurs 
whenever F() is a polynomial in 1/p without a constant term; in fact, the quadrature is exact 
for polynomials of degree not greater than 2n. G. SzEG6é [10] has shown that under quite general 
conditions a bret weaye type quadrature formula which converges for polynomials also converges 
for a much wider class of functions. Unfortunately his theorems do not seem to apply directly to 
the present case because the integral (1) involves a complex valued weight function which is not 
of bounded variation. 


1. H. S. Carstaw & J. C. JAEGER, Operational Methods in Applied Mathematics, 2nd edition, 
Oxford University Press, 1949, p. 75. 

2. H. E. Saurzer & R. Zucker, “Table of the Zeros and Weight Factors of the First Fifteen 
Laguerre Polynomials,’”’ Amer. Math. Soc., Bull., v. 55, 1949, p. 1004-1012. 

3. G. SzEG6, Orthogonal Polynomials, Amer. Math. Soc., Colloquium Pub., v. 23, 1939, p. 46-47. 

4. H. L. Kratt & O. Frinx, “A New Class of Orthogonal Polynomials: The Bessel Poly- 
nomials,’”” Amer. Math. Soc., Trans., v. 65, 1, 1949, p. 100-115. 

5. G. Dortscu, Theorie und Anwendung der Laplace-Transformation, Springer, Berlin, 1937, 


6. The shift in notation from (n+1) to m in A;™ will cause no confusion after the A;’s have 
been computed and are ready for use in (6). 
7. It was called to the author’s attention by H. L. Krai that P,(x) = (—1)*y,(x, 1, —1) 
where y,(x, a, b) are “‘generalized Bessel polynomials” (see [4]). 
. G. SzEG6, op. cit., p. 41-42. 
9. Formula (14) holds for » = 2 if we define Po(x) = 1. 
10. G. SzEG6, op. cit., p. 341-342. 


On the Improvement of the Solutions to a Set of 
Simultaneous Linear Equations using the ILLIAC 


The basic method used for solving simultaneous linear equations on the Uni- 
versity of Illinois’ electronic digital computer, the ILLIAC, has already been 
described in detail by WHEELER and Nasu [1 ]. The routine currently in use on 
the ILLIAC, programmed by Wheeler [2], makes use of the method of elimina- 
tion to solve the set of m simultaneous linear equations 


n—1 


(1) LX ax; + ain = 0 4=0,1,2,---,n—1 
i<0 


in a manner very similar to that used by a human solving such a system. 
In brief, the procedure used is as follows: 
a) The augmented matrix 


4=0,1,2,---,w—1 


2 aij 
2) ' j=0,1,2,-+-+,m 
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is scaled (by reducing each row by an appropriate factor) and punched as decimals 
with up to 12 places, row by row, on the input tape; each row is followed by 
an identifying symbol. 

b) The computer reads this matrix a row at a time, converts the numbers 
to binary form, and eliminates an additional variable from each additional row 
by using the previous rows which are stored in the memory of the machine. When 
this process is completed the matrix (2) is stored in ILLIAC as the triangular 
matrix 


4=0,1,2,---,n—1 


3 Cy 
(3) i j =0,1,2,---,m 
with ¢ = Ofor0 <j <i<n-—1. 

c) The solutions are obtained by a process of back substitution. Due to the 
fact that ILLIAC can represent numbers only in the range —1 to 1, provision 
must be made for automatic scaling of the results. If bars indicate the ILLIAC 
representation of a number (— 1 < # < 1), and if Z, = 10-” represents a scale 
factor, then the scaled solutions are formed during the back substitution by 
successively applying 
(4) ges SS jan —i,a—3%---4,0. 


fit Csi 
The true solutions are given in terms of the scaled solutions by 


z; z; ; 
(5) «=> * ios t= 0,1,2,---,#—1. 





The program first attempts to carry out the back substitution with p = 1. If at 
any stage the sum in (4) exceeds scale, p is increased by unity and the process 
is repeated. After the back substitution has been successfully completed the 
solutions #; are printed in columnar form followed by #, which serves to locate 
the decimal point in the column. Thus the machine actually operates on the set 
of equations 
n—-1 
(6) dX a8; + Ciné, = 0 4=0,1,2,---,n—1 
i=0 
whose solutions can be forced to remain within scale by proper adjustment of Z,. 
This method is extremely straightforward and is in fact quite similar to that 
which would be used to solve such a set of equations ‘“‘by hand”’. The chief utility 
of such a process applied to an automatic computing machine is that it minimizes 
the space necessary to store the matrix of coefficients in the equations. Any 
method that demands storage of the entire matrix a;; requires n(n + 1) ~ n?® 
memory locations. This elimination method uses the input tape as auxiliary stor- 
age and forms the reduced triangular matrix c,; as this tape is read in and thus 
requires only $n(n + 3) ~ $n*® memory locations. 
If the set of equations being solved is poorly conditioned in having a nearly 
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singular determinant, 
(7) Det (a;;) small 7,7 = 0,1,2,---,m —1, 


the above outlined process can yield solutions which are in error compared to the 
true solutions. Such errors can be aggravated by rounding-off in the many 
arithmetical steps of the process and by the truncation errors due to the finite 
register length of 39 binary digits. In order to obtain an improved set of solutions 
in such circumstances an iterative application of the method of residues has been 
programmed [3] for the ILLIAC. If x’; indicate the solutions to the equations 
(1) obtained by the ILLIAC and if corrections ¢; are defined by 


(8) xe x's + 
where the x; are the actual solutions of (1), then the corrections can be found from 
n—1 
(9) L ayej+6;=0 += 0,1,2,---,n —1 
j=0 
where 6; are the residues defined by 
n—1 
(10) 8, Do age’; + Gin t#=0,1,2,---,n—1. 
i—0 


Just as the original solutions were in error due to theffinite register length and 
the rounding-off errors in the arithmetic processes so are the e¢; in error. It could 
not be expected that a single application of equations (9) to find the corrections 
e; would yield the correct solutions to (1) by the use of (8). 

Hence, it has proved useful to apply the method of residues several times and 
also artificially to increase the numerical size of a residue contained in the machine 
by multiplying the residue by a factor B = 2* in order to gain numerical accuracy 
in the solution of (9). Let superscripts zero designate the solutions obtained by 
the conventional solution of the equation (1) as outlined in the above equations 
(3)—(5). Then if x; represent the solutions at the k** stage of the iterative process, 
the residues at this stage are defined by: 


n—1 
(11) 6x, = PD ans + ainx. + =0,1,2,---,n —1, 
<0 
where ILLIAC representations of the numbers are implied even though the bars 
have been dropped. The scale factor x, as determined by the first solution of 
equations (6) is maintained throughout the process. The corrections «;“ at the 
kt» stage are determined by dividing the result of (11) by x, to form 6; and 
solving the equations: P 


n—1 
(12) Lr age + B® =0 + =0,1,2,---,m —1. 
i=0 
The sum in (11) is formed by double precision continued multiplication which is a 


facility readily available on the ILLIAC. Hence the 6; are capable of accepting 
a multiplication by the scaling up factor B and still retain a register of significance. 
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The equations (12) are identical in form to equations (1) and can be solved in the 
same fashion. The (k + 1)** approximation to the solutions are: 


Xn 
x_*) 





(13) xD = x) 4+ Bre 4=0,1,2,---,n—1 


where x is the automatic scale factor supplied by ILLIAC in the solution of 
(12) which guarantees that each iterate is expressed in the same scale. 

To use Library Routine 100 the scaled coefficients of the augmented matrix 
are placed upon the input tape just as described in the use of Library Routine 51. 
To produce the x; the internal operation of the routine is identical to steps 
(a)—(c) outlined above. On each subsequent reading of the input tape the routine 
computes the residues (11) as each row of the input tape is read, solves the equa- 
tions (12), forms the improved solutions (13), and prints them. At the k stage 
the factor B is set to be 


(14) B= 2 %=1,2,---. 


If at any stage, multiplication by such a factor would cause the expression on the 
left hand side of equation (11) to exceed scale, the machine stops, replaces B 
by B’ = B2~ and puts itself into such a condition that the reading of the coeffi- 
cient tape for the kt stage can be repeated. If the equations are so badly condi- 
tioned that equation (11) will accept multiplication by no B > 1 without exceed- 
ing scale, this program is to no avail; the machine stops after producing a charac- 
teristic indicating signal. If the problem has been removed from the machine 
after the k* approximations to the solutions have been obtained it is unnecessary 
to repeat the process from the beginning if at a later time further iterations are 
desired. Provisions have been made so that the programnr, after being placed in 
ILLIAC, will first read in the x,“ from some previous iteration and store them 
in such a way that the condition of the machine is the same as if the iterative 
process had been continuing and had reached the k** stage. Subsequent iterations 
can hence be resumed from this stage. 

This method of annihilation of residues was chosen and used in ILLIAC to 
minimize storage requirements. By making the arithmetic procedures during 
subsequent iterations as nearly like those used on the first solution, the length of 
the operating routine stored in the memory has been reduced since most of it 
serves double duty. Again, only a triangular matrix need be stored at each itera- 
tion stage. This economy of memory space has however been achieved by the 
sacrifice of some speed since the tape of coefficients must be read on each iteration 
and the triangular matrix must be reformed at each iteration. The total tape 
reading time is at most yy5kn* sec. where k is the number of iterations necessary ; 
“the total computation time is s}ykn* sec., which is essentially dominated by the 
n® multiplications necessary to form the reduced triangular matrix. This program 
was designed to use only high speed storage so that even with the increased effi- 
ciency of storage yielded by this procedure 7 is still limited to moderate values so 
that the above times are not excessive in spite of the k repetitions. 

This method would be undesirable if enough storage were available in order 
to allow m to become large (order of 100, say) especially if that storage were rela- 
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tively slow (for example, a drum). In such a case the entire matrix should be read 
and stored or processed but once (reading time proportional to m* not kn*) and 
some iteration scheme requiring only n? multiplications per iteration used (compu- 
tation time proportional to kn? not kn*). This gain in computation speed is not 
quite so favorable as it may appear due to the greatly increased access time to 
the slower memory. 

The ILLIAC’s memory is sufficient to accommodate a system of 39 equations 
when used with Routine 51. The additional length of Routine 100 restricts to 
37 the number of equations that it can handle. With 37 equations the operation 
time of Routine 100 is about 4 minutes per iteration. For fewer equations the 
operation time is very roughly proportional to the square of the number of 
equations. 

The utility of this residue annihilating routine may be illustrated by an 
example comparing the results obtained with Routines 51 and 100 for sets of 
equations purposely constructed to be poorly condition to various degrees. 

The set of equations: 


(15) (1 —10-%) x1 —x2e+tx3—x,+3=0 
Xi — Xe +x3—x%+2=0 
X%1— Xe +1=0 
x1 —-x%+3=0 


has been constructed to have 


(16) | Det (a;;) 1 = 10-8 
1 = 102 
x2 = 10°+1 
x3 = 10%°+2 
x, = 10% + 3. 


These equations were solved to eleven decimal digits for 1 < gq < 9 by both 
Routines 51 and 100 as indicated in parentheses with the results: 





qt 
x; (51) x; (100) x: (100) 
10.000 000 05 10.000 000 004 10.000 000 000 
11.000 000 02 11.000 000 000 11.000 000 000 
12.000 000 06 12.000 000 004 12.000 000 000 
13.000 000 04 13.000 000 002 13.000 000 900 
qg=2 
2x;(51) x; (100) x; (100) 
100.000 000 20 100.000 000 20 100.000 000 00 
101.000 000 16 101.000 000 16 101.000 000 00 
102.000 000 20 102.000 000 20 102.000 000 00 
103.000 000 18 103.000 000 18 103.000 000 00 
q=3 
(51) x; (100) x; (100) 
1 000.000 000 1 1 000.000 000 2 0 999.999 981 8 
1 000.999 999 8 1 000.999 999 8 1 000.999 981 8 
1 002.000 000 2 1 002.000 000 2 1 001.999 981 7 
1 003.000 000 0 1 003.000 000 0 1 002.999 981 8 











2;(51) 


10 000.001 
10 001.001 
10 002.001 
10 003.001 


2x;(51) 
100 000.181 
100 001.181 
100 002.181 
100 003.181 


x: (51) 
1 000 018.190 
1 000 019.189 
1 000 020.190 
1 000 021.190 


x;(51) 


10 001 819.176 
10 001 820.172 
10 001 821.176 
10 001 822.174 


(51) 


100 182 149.38 
100 182 150.35 
100 182 151.39 
100 182 152.36 


820 
817 
821 
819 


9 


2 
9 
3 
1 


1 


i) 


q=4 
x; (100) 
10 000.001 
10 001.001 


10 002.001 
10 003.001 


qg=5 

x; (100) 
100 000.181 
100 001.181 


100 002.181 
100 003.181 


q=6 

x; (100) 
1 000 018.190 
1 000 019.189 
1 000 020.190 
1 000 021.190 


qgq=7 

x; (100) 
10 001 819.176 
10 001 820.173 


10 001 821.176 
10 001 822.175 


q=8 

x5 (100) 

100 182 149.38 
100 182 150.35 


100 182 151.38 
100 182 152.37 


xi (100) 


100 000 000.60 
100 000 001.61 
100 000 002.61 
100 000 003.61 


q=9 
x: (100) 


1 037 735 849.4 
1 037 735 850.1 
1 037 735 851.9 
1 037 735 852.1 


821 
817 
821 
819 


92 
88 


90 


meno w 


0 
0 
0 
0 





x (100) 


10 000.000 
10 001.000 
10 001.999 
10 002,999 


x; (100) 


10 000.000 
10 001.000 
10 002.000 
10 003.000 


xi (100) 


100 000.000 
100 001.000 
100 002.000 
100 003.000 


x4 (100) 


0 999 999.999 
1 000 000.999 
1 000 001.999 
1 000 002.999 


xi (100) 
000 000.000 


x; (100) 


09 999 999.669 
10 000 000.669 
10 000 001.669 
10 000 002.669 


x; (100) 


10 000 000.000 
10 000 001.000 
10 000 002.000 
10 000 003.000 


xs (100) 


099 999 668.22 
099 999 669.22 
099 999 670.22 
099 999 671.22 


x;® (100) 


100 000 000.00 
100 000 001.00 
100 000 002.00 
100 000 003.00 


x: (100) 


999 301 188.0 
999 301 189.0 
999 301 190.0 
999 301 191.0 
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999 


$3ss 333s 


ass 


oooo 
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x; (100) x; (100) 
1 000 026 370.3 0 999 999 511.6 
1 000 026 371.2 0 999 999 512.6 
1 000 026 372.3 0 999 999 513.6 
1 000 026 373.2 0 999 999 514.6 
x; (100) x; (100) 
1 000 000 018.4 0 999 999 999.3 
1 000 000 019.4 1 000 000 000.3 
1 000 000 020.4 1 000 000 001.3 
1 000 000 021.4 1 000 000 002.3 
x; (100) 
1 000 000 000.0 
1 000 000 001.0 
1 000 000 002.0 
1 000 000 003.0 


Several remarks about these results are in order. It is not to be expected that 
the zeroth order solutions of Routine 100 should agree with the solutions of 
Routine 51 since the conversion routine which reads the matrix into ILLIAC is 
slightly different in the two programs. The one in Routine 51 is accurate to 
+2-*, that in Routine 100 to +2-®. It is to be noted that as the equations 
become more poorly conditioned a greater number of iterations are necessary 
before Routine 100 converges to the correct answer. 

The deterioration of the Routine 100 solution for g = 3 upon iteration and 
the already excellent Routine 51 solution for g = 9 demand especial attention. 
Both of these effects are traceable to the rounding off of the scaling factors to 39 
binary digits in the machine. For the case q, a scaling factor 10-*~ is used by 
both Routine 51 and Routine 100 for the first iteration in accordance with equa- 
tion (5). In addition Routine 100 uses the scaling factor again on each subsequent 
iteration to form the residues and to reduce the resultant corrections to the 
uniform scale of the zeroth order solutions. The machine representation of the 
scaling factor was examined for each of the above cases. For q = 3, the scaling 
factor 10-* must be rounded down by six parts in the 11th sexadecimal digit 
when the number is rounded to ILLIAC’s register size of 10 sexadecimal digits. 
This is the largest rounding-off error encountered in the sequence. This accounts 
for the deterioration of Routine 100’s first iterate since the correction is multiplied 
by this scaling factor before being added to the zeroth order solution. For g = 9, 
the scaling factor 10- when represented to 10 sexadecimal digits is also exact 
to 11 sexadecimal digits, thus accounting for the already excellent Routine 51 
solution. Such effects would argue for the use of binary scaling factors during the 
carrying out of the back substitution. However, if the answers were presented 
with such factors they would not be immediately useful; if the answers were 
reconverted with the appropriate decimal scaling factors before presentation not 
only would additional memory space be required for such reconversion but the 
imprecision of these decimal scale factors would be introduced. 

Recent changes in the ILLIAC have resulted in an increased efficiency in the 
order alteration portions of a program. Various routines of the Program Library 
are being rewritten to utilize this increased efficiency. When Routine 100 is 
rewritten, sufficient memory space will be gained so that extra orders can be 
added to the program which will allow the residues to be computed more accu- 
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rately. This should speed the convergence and should eliminate most of the effect 
of imprecision in the scaling factor. 


JaMEs N. SNYDER 
University of Illinois 


Urbana, Illinois 
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The Use of Iterative Methods for finding the 
Latent Roots and Vectors of Matrices 


In a recent note in MTAC E. Bopewic [1] presented what he claimed was 
“a practical refutation of the iteration method for the algebraic eigenproblem.” 
In my opinion this note gave an entirely misleading impression of the value of the 
iterative method. Moreover an example was chosen as the basis of this refutation 
which so far from serving the purpose for which it was used, is in fact quite well 
suited to the iterative method provided it is used in a flexible manner. The 
iterative method, supplemented by a number of simple devices for accelerating 
convergence, has been used very effectively on the Pilot ACE to find the latent 
roots and vectors of a very large number of matrices, symmetric and unsymmetric, 
real and complex, up to orders as high as 60. Very high accuracy has been achieved 
even when many or all of the latent vectors of a matrix have been wanted. The 
details of the method used have been described in a recent-paper [2] by WILKIN- 
SON, but since that paper was written a magnetic drum store has been added to 
the Pilot ACE and the speed of iteration has thereby been considerably increased 
particularly for the larger matrices. The addition of the drum has also led to 
modifications in the details of the programme which have made it much more 
satisfactory to use. For this reason a general description of the programme is 
given below. The note includes an assessment of the value of iteration and con- 
cludes with the results achieved with it on Bodewig’s example. 

An understanding of the iterative programme will be aided by a description 
of the one or two facilities provided on the Pilot ACE, which are employed therein. 
There is a register, called the input register, which stores one of the standard 
words of 32 binary digits, into which a number may be inserted manually by 
means of 32 keys. The number thus inserted is displayed on a set of 32 lights and 
this number may be changed at any time during computation by the operator. 
The input register is addressable in the same way as all the other storage registers 
and the machine has access in 32 microseconds to the number held there, but it 
cannot send a word /o the input register. There is a second register, the output 
register, the contents of which are also displayed on a set of 32 lights. The machine 
may send a number to the output register in 32 microseconds, but it cannot read 
the number stored in it. The machine is also equipped with a monitoring device 
on which are displayed the contents of 32 consecutive storage registers. The 
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monitor may be made to display different groups of 32 words by the appropriate 
setting of a number of keys. The 32 words displayed on the monitor appear on 
32 rows in the form shown in fig. 1 where the first three lines only are shown with 
words of 10 binary digits instead of 32 for convenience. 
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The matrix A of which the latent roots are required is stored on the drum, 
and the present drum can deal with matrices up to order 60. Information may be 
transferred from the drum to the high speed store in blocks of 32 words at a time, 
each such transfer taking place in a time which, in this programme, is only 10% 
of the time taken to use that block of information. The machine starts with an 
arbitrary vector yo (the programme is arranged to read in a good initial guess 
if one is known) and from it forms two sequences of vectors y, and z, defined by 
the relations 


z, = (A — pl)y, 


z,/ (element of z, of maximum modulus) 


Ve+1 


where ? is the number set up on the input register. The normalisation of z, at 
each stage is carried out in order to keep numbers within range and it has the 
effect that the largest element in each y, is unity. Since p is read in at the be- 
ginning of each iteration only, it may be changed at any time during an iteration 
without affecting the value which is used for the remainder of that iteration. 
The vectors y, will tend to the latent vector, x, corresponding to that value of A 
(assumed real) which is either the algebraically largest latent root or the alge- 
braically smallest, according to the value of ». The value of » may also be chosen 
to improve the rate of convergence. Although is set up as a binary number this 
is not particularly inconvenient because there is no need to choose » with any 
great precision, so that most of the less significant digits are left as zero. The 
progress of the convergence may be studied by observing the successive approxi- 
mations to \ which appear on the output register. Since the time taken to send 
a number to the output register is 32 microseconds this display does not add 
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appreciably to the time of an iteration, as would printing successive values of X. 
As y, tends to the latent vector, the largest element of z, tends to (A — ») and 
the approximation to \ which is sent to the output register at the end of each itera- 
tion is found by adding # to this largest element. It is well known that for a sym- 
metric matrix the RAYLEIGH value given by 


‘ 
— Pe 
Ves 
has, in general, twice as many correct figures as the vector, y,, from which it is 
calculated, but no attempt is made to produce such a succession of \’s, the con- 
vergence of which would precede that of the vectors. This would be undesirable 
in that the convergence of the \’s would not then give a true impression of the 
convergence of the vectors. The sequence of \’s chosen reflects fairly accurately 
the convergence of the vectors, but experience has shown that there is a distinct 
tendency for the \’s to converge before the vectors. The behaviour of the vector 
may be studied by observing the monitor which is set so that it displays the 
components of the vector z, or 32 of the components if the matrix is of order 
greater than 32. The provision of this display does not add to the computing time 
and great use of it is made in connection with two devices for accelerating 
convergence. 

At the beginning of the iteration process, all the digits of all the components 
of z, will change with each iteration, but after a few iterations the more significant 
digits of z, will become constant and this is easily observed by the operator. Since 
the values of set up on the input register will always contain zeros in the lower 
digital positions, these positions are effectively free and may be used as further 
controls on the programme. The two least significant digits of the input register, 
which are called Pi and P2, are in fact used in the following manner. If at any 
time the P1 digit is set up on the input register then at the end of the current 
iteration, the machine takes the last three iterated vectors 


Ye-2, Yer. and y, 
and from them forms the vector Y, the i** component, Y*, of which is given in 
terms of the 7** components of the three y vectors by the relation 
yi= y' oy, a (y*,-1)? : 
y',2 — 2y*-1 + yi, 





This vector Y is referred to as the AITKEN vector derived from ¥,~2, ys—1 and y,. 
In general, as described by WILKINSON [2], the vector Y will be a much improved 
approximation to the latent vector. In order to be able to form the Aitken vector 
at any time, it is only necessary to store one extra vector. Suppose the successive 
vectors y, are stored in the set of registers A and the vectors 2, in the set of 
registers B. Then on completing an iteration, if an application of the Aitken 
process is not required, y, is transferred from A to a set of registers C and gz, is 
normalised to y,,;, which is then stored in A. If at the end of the next iteration 
an application of the Aitken process has been requested by the operator then the 
transfer from A to C does not take place and when z,,:, which is in B, is normalised 
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tO ¥.42 this is written over 2,41 instead of in A. Thus we have y, in C, y¥.4; in A 
and ¥,42 in B and this is sufficient for an application of the Aitken process. The 
Aitken process may be used as often or as little as the operator demands, though 
there must be at least three iterations between each application. 

The P2 digit is used to modify the operation of the programme as follows. 
If during an iteration a P2 digit is set up on the input register, then the machine 
stores away the current y, on the drum for use when the subdominant latent 
vector is being found. Suppose the initial vector yo is expressed in terms of the 
latent vectors x1, X2, ---, X, of A by the relation 


Yo = at + ate + +++ + ankn. 
Then, after s iterations, y, will be parallel to the vector 
ayAs'X1 + aede"X2 + +++ + GnAn'Xn- 


and after some iterations we will have in general 
@1A1° > aed2* > asr3° > °** D> andr’. 


This means that after a number of iterations the main impurity in y, will consist 
of the x, component and y, will be approximately 


x1 + €Xe. 


When the iteration for the dominant vector x; is complete a matrix of order one 
less is formed by a process of root removal and it is shown by WILKINSON [2] that 


(y. -_ %1) 


where y, is the vector which was stored on the drum, will, in general, be a very 
good initial approximation to the latent vector of the reduced matrix. The P2 
digit will be used once only during the process of iterating for one vector, and 
usually the time that is chosen is when about half of the binary digits in the 
successive yy, are constant. Iteration for each latent vector is continued until all 
the digits of y, are constant or the last few digits are going through a repeated 
cycle of values. This is most important because if we are to use successive root 
removals, the latent vectors must be as accurate as possible. The decision when 
iteration for a latent vector is complete is made by the operator and the pro- 
gramme is designed to accept a signal advising it when this is so. The process of 
reducing the matrix [2] then takes place and the machine begins to iterate with 
the reduced matrix for the next vector, determining its initial guess as de- 
scribed above. 

The last two techniques described are invaluable and are sometimes so success- 
ful that, starting from the good initial guess, 3 iterations plus one application of 
the Aitken process are sufficient for complete convergence. 

One or two comments on the role played by the operator might be made at 
this point. In the first place it should be stressed that at no time does the success 
of the programme depend critically on the operator’s performance. A failure to 
produce the best value of », to apply the Aitken process or to store away a vector 
at the optimum time merely means that one does not gain as much as one might 
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have done. It is undoubtedly true that an operator becomes rather more proficient 
with the programme after a little experience with using it but the operators em- 
ployed on the Pilot ACE are in no sense of the word skilled mathematicians, and 
I think it is generally agreed that the programme is much more entertaining than 
any other used on ACE. 

In my opinion the “algebraic eigenproblem”’ is not a problem for which there 
is a universal solution. There are, in practice, many different requirements. The 
matrix may be symmetric or unsymmetric and in the latter case, the roots, real 
or complex. All the roots may be wanted, several dominant roots, the few smallest 
roots or the roots, if any, in a given range. Again we may or may not want the 
vectors. Since all machines have stores of finite size often divided up into high 
speed and auxiliary sections, storage considerations often have a vitally important 
part to play. 

All possible combinations of the above requirements have been met in using 
the Pilot ACE and the iterative process has proved to be about the most useful 
weapon at our disposal but it is certainly not the only method we use. For example 
if all the roots, but not the vectors, of a symmetric matrix are wanted the JACoBI 
process [3] is undoubtedly better than iteration unless the separation of roots is 
quite exceptionally good. If a few isolated roots of a symmetric matrix, but not 
the vectors, are wanted, then the method described by W. Givens [4] is almost 
certainly the most effective of known methods. The iteration method has been 
used mainly on unsymmetric matrices and most commonly where several of the 
dominant roots and the vectors have been wanted but it would be used equally 
on symmetric matrices where the requirement was the same. 

As regards the separation of successive roots, Bodewig’s statement that a 
separation of four to one is necessary is very wide of the mark. For some of the 
larger matrices, of orders 50 to 60, which we have received from aircraft firms, 
separations of one part in ten have been quite common and much worse separa- 
tions have been dealt with satisfactorily using the above accelerating techniques. 
Among smaller matrices quite bad separations have been treated very success- 
fully. Below are given the roots of three matrices A, B and C for which separation 
is bad. These three examples are chosen from a group of 30 unsymmetric ma- 
trices of orders 10 to 15 all of which gave separations comparable with those in 
the examples. 


5.5695580 


Roots of Matrix A Matrix B Matrix C 
0.2144582 0.2155673 0.1562830 
0.4644432 0.4826537 0.3571287 
0.6813217 0.7523862 0.9520832 
0.9824573 1.0217583 1.2017653 
1.2457623 1.5162837 1.8635762 
1.7428765 2.0256132 1.8801257 
2.1281317 2.6578326 2.4613258 
2.6315722 4.1243587 2.9980751 
3.2458126 4.5625830 3.5287125 
5.4916793 5.5612571 4.0165324 
5.5491920 5.6013757 4.5238175 
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In matrix A, for example, the three dominant roots are separated by one part in 
270 and one part in 100 respectively. It was immediately obvious that convergence 
for the dominant root was very slow so that an attempt was made to find the 
roots in the reverse order by an appropriate choice of p. A value of » was chosen 
which was appreciably greater than half the value of the initial approximations 
which had been sent to the output register when iterating with = 0. The values 
of » which were used were not kept but it is probable that the root 0.2144582, 
which was found first, was found by iterating with (A — pI) where p was be- 
tween 3 and 3.5. The separation between the two dominant roots of (A — 3J), 
for example, is about one part in ten and this with the accelerating process is quite 
adequate. The three largest roots were found last and by this time the reduced 
matrix was of order 3, so that iterations were taking place at the rate of 30 per 
second. Choosing p to be near the approximations to \ which appeared in the 
output register also gave quite good separation so that the largest three roots 
were found very quickly. For matrix C the two roots 1.8635762 and 1.8801257 
were found last. 

It is popularly believed that the eigenvectors obtained after several root 
removals tend to be fairly inaccurate. In this connection our experience on the 
Pilot ACE with several hundred matrices (literally several thousand if matrices 
of order less than 10 are included!) on which the iteration method has been used 
is most interesting. If we define the residuals corresponding to a latent root \ 
and a vector x as the components of 


Ax — «x 


then in no example has any residual been greater than 20 X 2-** (i.e., approxi- 
mately 37 X 10-*) where it is assumed that the elements of A are of order unity. 
Moreover more than 90% of all residuals are in fact below 5 X 2~-** and a set of 
latent vectors with all its residuals consisting of either zeros, +2-*°, +2-*8 is 
far from uncommon. The main reasons for the accuracy are that iteration for 
each vector is taken to the limit and scalar products are accumulated with all 
the digits produced by each multiplication (N.B., the elements of the matrix A 
and the vectors y, are single length numbers). No other method which we have 
used has produced such consistently accurate vectors. Processes in which the 
latent roots are found from the characteristic equation have, in our experience 
proved to be deceptively satisfactory when dealing with small matrices. For larger 
order matrices it very easily happens that such precautions are necessary in 
forming the characteristic equation that what appears to have been a very fast 
and effective method for small matrices, proves unexpectedly tiresome for those 
of higher order. Determining the latent vectors accurately, particularly in the 
case of unsymmetric matrices is even more hazardous, and it is therefore very 
misleading to try to assess the relative efficiency of methods by counting the 
number of multiplications they require. Bodewig remarks that it is possible to 
find the latent roots and vectors to any accuracy from the characteristic equation 
but it should be remembered that the accuracy is in fact limited by that of the 
characteristic equation which has been calculated and that when roots are close 
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together a given accuracy in the characteristic equation does not produce the 
same accuracy in the roots. 

It is interesting to see how the Pilot ACE programme for iteration dealt with 
the example given by Bodewig. The coefficients of the matrix were stored in the 
machine to a precision of 29 binary figures (or 8.7 decimals) for a reason which is 
associated with the details of the programme. Iteration was started with p = 0 
and it was immediately obvious that convergence, which one would expect to be 
almost instantaneous for a matrix of order 4, was in fact proceeding comparatively 
slowly. As a trial value, to speed convergence, » = 2 was set up on the input 
register. Convergence then took place almost immediately so that no thought 
was given to finding a better value of ». All four roots were found in a few seconds 
iteration time and the roots and vectors, after normalising to agree with Bodewig’s 
method, were as given below. 


2 1 3 4 
i =—3 1 5 
Pm 3 1 6 -—2 
4 Ss =—2) —1 


Aa = —8.02857835 Az = +7.93290471 As = +5.66886437 Ay = —1.57319073 


X1 


X2 X3 X4 
+1.00000000 +1.00000000 +1.00000000 +1.00000000 
+2.50146029 +0.37781815 +0.95700150 —0.90709211 
—0.75773064 +1.38662122 — 1.42046822 —0.37759122 
— 2.56421169 +0.34880573 +1.74331690 —0.38333124 


The vectors agree, apart from the end figure of one or two components, with 
those given by Bodewig, except for the final component of x, where the two re- 
sults differ by 0.00001886. It is easy to satisfy oneself that Bodewig’s result is in 
error, and it is probable that his results were copied incorrectly. The only peculi- 
arity displayed by Bodewig’s example was that its two dominant roots were almost 
equal and opposite in value and this is not a difficult situation. Indeed even if 
the roots had been exactly equal and opposite in value there would have been 
no difficulty whatever. If the two dominant roots had-had the values in the 
example but with the same sign the situation would have been just a little less 
satisfactory. Then the root —1.573 etc., would have been found first and the root 
5.669 second. The matrix obtained after two reductions would have been of order 
two, so that iteration would have been extremely fast and a choice of » anywhere 
near 8.0 would have given almost immediate convergence. Such a choice of p 
would have been suggested by the approximations to \ obtained on the output 
register. Besides being a poor choice on which to base a refutation of the iteration 
process, Bodewig’s example has other weaknesses which make it a rather poor 
one for estimating the relative efficiency of methods. In the first place it is of very 
low order so that it does not reveal the difficulties of preserving accuracy which 
arise with larger matrices and secondly, it has small integer coefficients. This 
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means, for example, that if we use a method in which we find the characteristic 
equation we will almost certainly obtain the latter exactly, which will not be true 
if the coefficients were numbers with several digits. 

J. H. WiLkinson 
National Physical Laboratory 
Teddington, Middlesex 
England 


This note is published with the permission of the Director of the National Physical Laboratory. 


E. Bopewie, “A cae refutation of the iteration method for the algebraic eigenproblem,” 
MTAG v. 8, 1954, p. 237-23 
Yee 4 WILKINSON, ‘ The ‘calculation of the latent roots and vectors of matrices on the Pilot 
Model a the ACE,” Cambridge. Phil. Soc., Proc., v. 50, 1954, p. 536-566. 
i. &, 3 Jacopr, ‘ ‘Uber ein leichtes Verfahren, die in der Theorie der Sdkularstérungen 
vaiheasaaniil Gleichungen numerisch aufzulésen,”’ J. reine angew. Math. 30, 1846, p. 51-94. 
4. W. Givens, “Numerical computation of the characteristic values of a rea ‘symmetric 
matrix,” Oak Ridge National Laboratory. ORNL1574. 
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The Specific Heat Function for a Two-Dimensional Continuum 
Numerical Values of 


C2 6 (7 #dét 2x 


C. eJoe—-1 e-1 





This function which appears in the theory of low-temperature specific heats 
of two-dimensional (layer) structures [1] was computed as follows: 


(1) For 0 < x < 2.0, the formula 


C2 x? x* xs x8 x0 
co Pike 24 - 720 24,192 + 864,000 31,933,440 + 





was used. The maximum error (using seven terms) is no greater than 0.5 X 10-5. 
(2) For 2.0 < x < 16.0, the formula 


C; _14.424684 » ~=(4 2 2 ) _ _ 2 
Gack x? Ox 2 ‘ nx + (nx)? + (nx)? e—1 


was used. The maximum error was approximately 2 X 10-*. 
The value C,, = 3R = 5.9616 cal mol deg was used throughout [2]. 


x C2/Cw Cz 
0.0 1.00000 5.9616 
0.1 0.99958 5.95911 
0.2 .99833 5.95167 
0.3 -99625 5.93925 
0.4 -99333 5.92186 
0.5 .98959 5.89955 
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C2/Co 


-98502 
.97961 
.97340 
-96632 
-95843 


95155 
.94276 
93336 
.92337 
-91284 


.90179 
-89026 
.87829 
.86592 
85318 


84024 
82677 
81327 
-79936 
-78526 


77111 
75689 
-74243 
-72802 
-71353 


.69901 
68454 
.67014 
65582 
64159 


.62746 
61347 
59961 
58594 
.57240 


-55910 
54597 
.53306 
52037 
50791 


49568 
48368 
47193 
46041 
44915 


C; 


5.87227 
5.84006 
5.80300 
5.76081 
5.71378 


5.67274 
5.62037 
5.56432 
5.50478 
5.44198 


5.37610 
5.30737 
5.23601 
5.16224 
5.08629 


5.00919 
4.92890 
4.84841 
4.76548 
4.68138 


4.59707 
4.51225 
4.42607 
4.34016 
4.25379 


4.16719 
4.08097 
3.99511 
3.90974 
3.82488 


3.74068 
3.65725 
3.57461 
3.49313 
3.41244 


3.33315 
3.25483 
3.17787 
3.10225 
3.02793 


2.95502 
2.88352 
2.81343 
2.74480 
2.67766 
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x C2/Co C; 
5.1 43814 2.61202 
5.2 42738 2.54789 
3.3 41687 2.48521 
5.4 .40661 2.42403 
5.5 .39659 2.36434 
5.6 .38683 2.30611 
§:7 .37731 2.24937 
5.8 .36803 2.19404 
5.9 .35900 2.14019 
6.0 .35019 2.08771 
6.1 .34163 2.03664 
6.2 .33328 1.98691 
6.3 .32517 1.93853 
6.4 .31728 1.89152 
6.5 .30961 1.84577 
6.6 .30215 1.80129 
6.7 .29490 1.75805 
6.8 .28785 1.71604 
6.9 -28100 1.67519 
7.0 .27434 1.63554 
7.1 .26788 1.59699 
7.2 .26160 1.55954 
7.3 .25550 1.52318 
7.4 .24959 1.48797 
7.5 .24382 1.45353 
7.6 .23823 1.42021 
7.7 .23280 1.38786 
7.8 .22753 1.35643 
7.9 .22241 1.32591 
8.0 .21744 1.29628 
8.1 .21258 1.26733 
8.2 .20792 1.23954 
8.3 .20337 1.21239 
8.4 .19894 1.18602 
8.5 .19465 1.16041 
8.6 .19048 1.13554 
8.7 .18642 1.11137 
8.8 .18248 1.08789 
8.9 .17866 1.06508 
9.0 .17494 1.04291 
9.1 .17132 1.02136 
9.2 .16781 1.00043 
9.3 -16440 -98009 
9.4 .16108 .96030 
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x C2/C. C; 
10.0 -14302 .85262 
10.5 -13006 .77535 
11.0 -11873 .70779 
21:5 .10897 .64965 
12.0 .09998 .59604 
12.5 .09220 .54965 
13.6 .08528 .50840 
13.5 .07910 47157 
14.0 .07357 43857 
14.5 .06859 .40890 
15.0 .06411 .38218 
p .06003 .35789 
16.0 .05634 .33589 

M. DANK 


S. W. BARBER 
Owens-Illinois Glass Co. 
Toledo, Ohio 


1. V. V. Tarasov, ‘“Teoriia Teploemkosti Tsepnykh I Sloistykh Struktur” [The Theory of 
Heat Capacity of Chain and Linear Structures], Zhurnal Fiz. Khimii, XXIV, 1950, p. 111-128. 

2. J. W. M. Dumonp & E. RicHarp CoHEN, “Least-squares adjusted values of ‘the atomic 
constants as of December, 1950,” Phys. Rev., v. 82, p. 555-556. 


A Continued Fraction for e 


Let A,/B, denote the convergents to the continued fraction 


ae an 


b+— . + eatte F + ..- 
It is well known [1] that the continued fraction 


abe 20 3b, AA sbobe5 
bibs + G2 ~ (bobs + as)bg + boa4 ~ (babs + as)be + dads 
AcA7b bg 
~ (beb7 + a7)bg + boas ~ 








bo + 





has convergents A on/Bon. 
From the continued fraction of Gauss, one can obtain [2] 





x 2f/1-2 ef23'' 2/25 2x2/2-5 
wr r 
hor in7 1 * eee tf + 1 
For this expansion, we have bo = 1, a; = x, Gon = — x[2(2m —1)], Gongs 


= x/[2(2m + 1)], and b, = 1 (mn = 1,2,3,---). Thus, we obtain 


— Arn@on41D2n—2don42 = x*/[4(4n? — 1)], 
and 


(DonDon+1 + Gon+1)D2n+2 + Don@on+2 = 1, (n = 2,3,:- *); 
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and so 


x x?/(4-3) x*/(4-15) 
i-2/2t 1 F 1 5 
x*/{4[4(» — 1)? — 1]} 
+ 1 + 


e=1+ 








A convenient equivalent form is 


ee eg 2/4 
trlittTri+ a e's We 4 Se 


These expansions converge quite rapidly forallx. Forexample, if — 1 < x < 1, 
we have e — Au/By < .000084, e — A;/Bs < .000000033, and e — Ac/Bs 
< .000000000081. 

N. Macon [3] 
Alabama Polytechnic Institute 
Auburn, Alabama 


This note was written while the author was a summer participant at the Oak Ridge National 
Laboratory, 1954. 


1. See, for example, OskKAR PERRON, Die Lehre von den Kettenbriichen, 2nd ed., Teubner, 
Leipzig, 1929, p. 201. 
2. PERRON, op. cit., p. 312. Mr. Garwick also notes that this is credited to Thiele, see, for ex- 


ample, N. E. N¢grlund, Vorlesungen siber Differenzenrechnung, Springer, Berlin, 1924, Chapter 15, 
p. 415-455, especially p. 454. 


3. The referee has noted that this work has been done independently by J. V. Garwick of the 
Norwegian Defense Research Establishment. 
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78(C, D, E].—Friepricu Léscu, Siebenstellige Tafeln der elementaren transzen- 

denten Funktionen. Springer, Berlin-Géttingen-Heidelberg, 1954, viii + 335 p. 

27 cm. DM 49.80. 

This volume is meant to supersede the 1926 table of K. Havasai [1], now 
out of print. The elementary transcendental functions are sin x, cos x, tan x, 
sinh x, cosh x, tanh x, In x, e*, e~*, arcsin x, arctan x, arcsinh x, arctanh x, together 
with yg, the angle corresponding to the arc x. In Table I (40 pages) these are 
given to 9D for x = 0(.0001)0.1, together with first differences (except for In x, 
and, of course, g). In Table II (277 pages) these are given to 7D for x = 0.1(.0005) 
3.15, with (double) first differences for all functions (but ¢) and for x = 3(.01) 
10(.1)20, without differences. For x > 1, arccothx is tabulated in place of 
arctanh x and arccosh x in place of arcsin x. There are supplementary tables of 
tanx to 7D for x = 1.5680(.0001)1.5730, arctanh x for x = 0.9980(.0001)1, 
arccoth x for x = 1(.0001)1.0020. 

In Table III, there are given sin x, cos x, In x, arcsinh x, arccosh x to 7D, 
e+* to 7S, for x = 0(1)100. Table IV gives $n to 12D for m = 0(1)100. Table V 
gives sin rx, cos $ax to 7D for x = 0(.001)0.5. Table VI gives exp (+422), 
sinh 44x, cosh 44x to 7D for x = 0(.01)2. Table VII gives exp (+2x/180), 
sinh (#x/180), cosh (rx/180) to 7D for x = 0(1)180. Table VIII gives ¢ in degrees 
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(12D) and in degrees, minutes and seconds (8D) for x = 0(.0001).001(.001) 
01 (.01).1(.1)1(1)10(10)100. Table IX gives x in radians to 12D for g = 0(1°)90°, 
for ¢ = 0(1’)1° and for g = 0(1”’)1’. Table X gives a collection of constants to 
15D, the exact values of the first twenty Bernoulli numbers and the first 15 
factorials; also (mm!) to 15S, logio ! to 15D and the binomial coefficients (”) 
for m = 1(1)15,7 = 1(1)[4n]. 

The volume is excellently printed on rather large pages; it is, however, ex- 
pensive. The Introduction gives worked examples, which cover typical cases. 
In general the error in linear interpolation (in Table I, II) is at most two units; in 
cases where it is larger, attention is called by use of italics for differences or by 
their omission. 

It is stated that considerable parts of the tables were computed especially for 
this volume ; for the remainder use was made of available tables and interpolation 
in such tables; the calculations were made to 9 or 10D and checked by differencing 
before rounding to 7 or 9D. The new material is mainly that concerned with the 
inverse hyperbolic functions and the inverse tangent; some of the hyperbolic 
material is available in the Harvard volume [2 ]. 

Two pages of Table I and two of Table II were checked by comparison with 
recognized sources and new calculations (made on SEAC): no discrepancies 
were noted. 

Users will have to be wary of noting the leading figures of values and differences 
at the head (and foot) of the columns and the changes which occur in them. 


l *f » 


1. K. Havasu, Sieben- und mehrstellige Tafeln der Kreis- und Hyperbel-funktionen und deren 
Produkte sowie der Gammafunktion. Springer, Berlin, 1926. 

2. HARVARD UNIVERSITY, Computation Laboratory, Annals, vol. 20, Tables of inverse hyper- 
bolic functions, 1949. 


mk 


79[D, I].—Ow_En R. Mock, A table of sines and cosines of 9256" 


13 typewritten 


leaves deposited in UMT FILe. 

These tables, constructed for Fourier analysis and synthesis with data from 
binary sources list sin and cos to 25D for k = 1(1)255. There is no introduction 
or description. The table is believed to have been computed on a Model 701 
International Business Machine. The author believes that all digits are correct. 


CoB T. 


80[F ].—Emma Leume_r, “On the cyclotomic numbers of order sixteen,’’ Canadian 
Jn. Math., v. 6, 1954, p. 449-454. 


Dickson proved that (i, 7)s, the number of solutions of 


girtt ok i= g®#ti(p) 


can be expressed for every i, j as a linear combination, with rational coefficients 
whose denominators are divisible only by factors of 64, of p, x, y, a, b, where 
p = x? + 4y’ = a? + 20’ and pis a prime of the form 8” + 1. The signs of x, y, a, b 
have to be chosen in a special way which was discussed in another paper by the 
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author “On the number of solutions of u* + D = w* (mod p),” Pac. J. Math., 
v. 5, 1955, p. 103-118. 

In order to decide whether analogous relations hold for (i, j)1s, with coeffi- 
cients whose denominators are divisible by factors of 256, computations were 
carried out on SWAC for a set of prime numbers of the form 32m + 1 and 32n + 17. 
These computations lead to the conclusion that Dickson’s result cannot be ex- 
tended to 16, at least not with the same sign rules which worked for (i, j)s. 

Some of the numerical results obtained for this investigation are published 
here as a table. It is possible to obtain 205 of the 256 values of (i, 7)1s from the 
51 values of (7, j)1s for which it < j, 7 —i > 1, 16 — j >i. Only these (3, f)i¢ 
are tabulated for the eight primes of the form 32m + 1: 


97, 193, 449, 641, 673, 769, 929, 1409 
and for the six primes of the form 32n + 17: 
241, 401, 433, 977, 1009, 1297. 


Also the values of x, y, a, 6 are given for each p, with the signs needed for the 
relations concerning (i, j)s. 
Printing error: p. 449, line 7 from top should read: x = a = 1(4). 
O. Taussky 


National Bureau of Standards 
Washington, D. C. 


81[F ].—Francis L. Misa, “‘A table of integral solutions of a? + &? + 2 = 7°,” 
Mathematics Teacher, v. 48, 1955, p. 251-255. 
Solutions for r = 3(2)207. See Review 82 in this issue. 
This paper lists 14 errata which the author has found in the earlier table by 
R. H. Bacon [1], which listed solutions for r = 3(2)99. 
am Ty A 


. Bacon, “Integral solutions of x* + »* + 2* = r*,”" School Science and Mathematics, 
Vv. 4, 1947, p. 155-164. 


82[F].—Francis L. Mixsa, A table of integral solutions of A* + B* + C* = R® 

for all odd R from R = 3 to R = 325, in two sets, R = 3 to 207 and an exten- 

sion from R = 207 to R = 325. 32 leaves, reproduced from typewritten 

manuscript. On deposit in UMT FILe. 

This table extends earlier tables by R. H. Bacon [1] and the author [2] to 
R = 325. Bacon’s table listed decompositions of odd numbers into the sum of 
three squares for numbers R = 3(2)99, and the paper by Mrxsa presented the 
list for R = 3(2)207. The calculation follows the lines of the two described; the 
present table has a clear introduction telling how it was calculated and listing a 
few pertinent identities. For R > 209 non-primitive sets are marked. 

©. 3.2% 


1. R. H. Bacon, “Integral solutions of x* + 7 + 2* = r*, School Science and Mathematics, 

Vv. @. S47, Pp. 155-164. 
” Francis L. Mrxsa, “‘Table of integral solutions of a? + b* + c = r*,” Mathematics Teacher, 

v. és, 1955, p. 251-255. 
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83[F].—Francis L. Mrxsa, Table of quadratic partitions x* + y? = N. Manu- 
script of 102 leaves (reproduced from typewritten masters) deposited in the 

UMT Fre. 

The table attempts a complete description of the possible decomposition of all 
odd numbers N = 1(2)50,000 into the sum of two squares. Some pertinent in- 
formation concerning this is contained in [1 ]. Rules for deriving similar decompo- 
sition of even numbers are given in the introduction. No details of the calculation 
are included. 
or eT: 


1. J. V. Uspensxy & M. A. HEASLET, Elementary Number Theory, McGraw-Hill, New York, 
1939, p. 337-341. 


84[F].—R. J. Porter, Irregular Negative Determinants of Exponent 3n with their 
critical classes. Part I, from — D = 1 to 50,000. 174 typewritten pages, 25.5 
X 10 cm., deposited in UMT Fine. 
This table is a supplement to Table of irregular negative determinants of ex- 
ponent 3n up to — 2 = 150,000, which were previously deposited in the UMT 
Fite [MTAC, v. 9, 1955, p. 26]. 


R. J. PORTER 
266 Pickering Road 


Hull, England 


85[F, I, J, K].—Francis L. Miksa, Stirling numbers of the second kind. Type- 
written manuscript (reproduced from typewritten masters) of 32 p., on deposit 
in UMT Fire. 

This table gives ,.S, for m = 1(1)50, a table of G, for m = 1(1)51, an excellent 
introduction, and a bibliography prepared by HENRY Frn_ayson. The number 
-S, is the Stirling number of the second kind defined by the expansion formula. 
It is also given as 

A’0" 
“a 
The number G, = >, -Sn. 

The author gives the recurrence relations used in computing and in checking 
the tables. The numbers are used in number theory, some aspects of analysis 
dealing with series and with probability or distribution problems involving ar- 
rangements of objects. 

Other published lists are reported by FLETCHER, MILLER and ROSENHEAD [11 ]. 

es ¢ ¥ 


1. FLETCHER, MILLER & ROSENHEAD, An Index of Mathematical Tables, Scientific Computing 
Service, Ltd., London, p. 72. 


86[H_].—Anpres Zavrotsky, Tablas para la Formacion de las Ecuaciones Cubicas, 
Publicaciones de la Direccion de Cultura de la Universidad de los Andes, 
No. 44, Merida, Venezuela, 1955, iv + 48 pp., 21.5 K 15.8 cm. 
This table gives the exact values of the coefficients of the cubic equation 
x? + px + g = 0asa function of two of its roots x; and xz (the third root x; being 
— x1 — x2). The functions —p and g are tabulated side by side for x, = 1(1)100 
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and x2 = x:(1)100, where — p = x;? + x2 + xix2 and q = x %2(x; + x2). There 
is a brief introduction containing some notes on the early history of the subject 
and some obvious comments on how to obtain » and g from the table when x; is 
either negative or non-integral, which is followed by some dozen references to 
other tables for solving cubics. 


HERBERT E. SALZER 
Diamond Ordnance Fuze Laboratories 


Washington, D. C. 


87(H, I].—Epvarp L. St1eFeL, Kernel polynomials in linear algebra and their 
numerical applications, multilithed from typescript, 1955, 52 p., 26.7 X 20.3 
cm., possibly available on loan from libraries or the author. 

High-speed digital computing machines make it possible to deal with matrices 
of much larger orders than would have been considered fifteen years ago. As a 
consequence, in the last decade many older matrix methods have been re- 
examined, and a number of new methods—mostly iterative—for treating linear 
problems have been devised by FRAME, HESTENES, KANTOROVICH, LANCzOs, 
PURCELL, SOURIAU, STIEFEL, D. YOUNG, and others. The practical computer faces 
an almost impossible task of remembering all these methods for solving linear 
problems, not to mention understanding them or being able to predict which is 
best for a given task. The reviewer's outline [1 ] of methods of solving linear sys- 
tems is terrifying by its length and lack of unity. 

But now Professor Stiefel has unified several of the known iterative methods 
for matrix probiems. The unifying principle is the theory of orthogonal poly- 
nomials over an arbitrary mass distribution on an interval of the real axis. The 
booklet reviewed here is a clear and careful exposition of this material. 

In the first section the author reviews the theory of orthogonal polynomials 
{P;(A)} over a mass distribution p(A)dA. He shows how one can generate them by 
orthogonalizing the sequence {AP;_:(A)}, a far more efficient process than or- 
thogonalizing the powers {\‘}, as Lanczos has stressed. Although the author 
doesn’t mention it here, this idea can be worked into a polynomial curve-fitting 
routine much more effective than the usual ones based on solving normal equations 
built from {d*}. 

To solve a linear system Ax = k, the author considers simple iterative rou- 


tines of the form x41 = x; + Axi, where Ax; = (1/q,)ri, 7; = k — Ax;. Then 
= R,(A)k, where R,(A) = Tl (1 — A/q) is a polynomial with R,(0) = 1. 
r=0 


The author poses the following problem of choosing a best strategy: Suppose the 
eigenvalues of A are known to lie in an interval [a,b]. Then it is desired to 
find that polynomial R,(A) of maximal degree m such that R,(0) = 1 and such 
that f/.° R,(A)*p(A)dA is minimized. Since one doesn’t know the eigenvalues ),, 
such a requirement is about all that one can make about the smallness of R,.(A,). 
The author solves this minimization problem exactly in terms of orthogonal 
polynomials, and shows various ways in which it may be used to construct 
algorithms which yield small r,. He points out, however, that the intermediate 
residuals r; may be quite large; i.e., these processes may be unstable. 

The author now shows how to modify the above algorithm so that each inter- 
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mediate r; also corresponds to a best strategy, thus restoring stability. The change is 
to put Ax; = (1/q:) (ri: + piAx:_1), where p;, g; are scalars involved in the three- 
term recurrence relation for {P;(A)}. He shows that various known algorithms 
correspond to special choices of p(A)—in particular, the conjugate gradient method 
of Lanczos, Hestenes, and Stiefel. 

In another section the same ideas are applied to the calculation of eigenvalues 
of a matrix A with real eigenvalues. It is shown how to bring out any intermediate 
eigenvalue by iterations with an appropriate kernel polynomial. As a special case 
the author shows us the intriguing new spectroscopic method of Lanczos. The 
reviewer feels that the ideas of this booklet will prove more useful in eigenvalue 
problems than in solving linear systems, because of the greater number of other 
ways of dealing with linear systems. 

The final section is devoted to an application of the quotient-difference al- 
gorithm of RUTISHAUSER and Stiefel to generating the orthogonal polynomials. 

This booklet is an exceedingly valuable contribution, and it is hoped that the 
material will be made available to every mathematician concerned with numerical 
methods in linear algebra. A person studying the booklet carefully will see— 
probably for the first time—the unity behind the gradient methods, RICHARDSON’s 
method, and the various conjugate gradient methods. And he will learn how to 
improve his techniques of curve fitting. 

The author uses the notation of the Dirac delta function for atomic mass 
distributions, and also to motivate the derivation of the kernel polynomials. 
This is most unusual for a mathematician, and is probably done in an effort to 
reach non-mathematical readers. While the reviewer is neutral on the matter, 
he finds that several of his colleagues resent it emphatically. The reviewer’s only 
criticism is that there is no direct reference to the work of L. F. Richardson, 
who in his 1910 paper [2] was considering similar questions of strategy. The 
reviewer also wonders whether the author couldn’t as well assume at many places 
in the exposition that A has complex eigenvalues, and thus broaden the treatment 
substantially. After all, a limitation to real eigenvalues is almost a limitation to 
symmetric matrices. 


GEORGE E. FORSYTHE 
University of California 
Los Angeles, California 


4.:G, E. FORSYTHE, “Tentative classification of methods and bibliography on solving systems 
of linear equations,’ " National Bureau of Standards Applied Mathematics Series 29, 1953, Superin- 
ae ae 4 of Documents, Washington, D. C., p. 1-28. 

L. F. RicHarpson, “The approximate arithmetical solution by finite differences of physical 
olin we differential uations, with an application to the stresses in a masonry dam,” 
Roy. Soc. London, Phil. Trans. (A A), v. 210, 1910, p. 307-357. 

Editorial note. These were notes of four lectures prepared under contract between the National 
Bureau of Standards and the American University under sponsorship of the Office of Naval 
Research. They will probably be published as part of a volume to appear within a few months 
as part of the NBS Applied Mathematics Series. 


88(H, K, P, S, Z].—Transactions of the Symposium on Computing, Mechanics, 
Statistics, and Partial Difference Equations held at the University of Chicago, 
April 29-30, 1954. Vol. II. Editorial Committee: F. E. Grusprs, F. J. Murray, 
& J. J. Stoker, Interscience Publishers, Inc., 1955, 216 p., $5.00. 

This volume contains eleven papers based on invited addresses given by the 
separate authors at the second symposium on applied mathematics sponsored 
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jointly by the American Mathematical Society and the Office of Ordnance 
Research, U. S. Army. This symposium was held on April 29 and 30, 1954 at 
the University of Chicago. A great variety of different fields of applied mathe- 
matics and mathematics is displayed in this book, ranging over operational 
research, statistics, and partial differential equations. 

P. M. Morse discusses Operations Research as an effort to discover regularities 
in the behavior of teams of men and equipment as they go about doing assigned 
jobs and to link these regularities with other knowledge so that the phenomena 
occurring can be modified and controlled. Examples of the application of opera- 
tions research to various problems are given. 

J. NEYMAN gives a consistent presentation of his fundamental contributions 
to the theory of inductive inference and answers objections raised by CARNAP. 
Neyman propounds the view that prior to any discussion based on a theory of 
inductive inference a theorizing step must be taken to create a mathematical 
model of the phenomenon being studied. He points out that the model may differ 
from the theory of probability based on the axioms of Kolmogoroff. The general 
discussion of the first part of the paper is applied to the problem of homogeneity 
of neutral V-particles. 

H. O. HARTLEY reports progress on two topics in his paper entitled, ‘‘“Some 
Recent Developments in Analysis of Variance.’”’ These topics are: (1) multiple 
decisions and comparisons and (2) short cut procedures based on substitute 
measures of dispersion. 

The ‘“Two Unsolved Problems of Statistical Mechanics” explained by J. E. 
MAYER are: 


(1) The evaluation of the pressure P of a system with N molecules as a func- 
tion of the activity Z and the temperature T as a sum over integrals of functions 
of the potential energy over configurations in 3 N dimensional space. The equiva- 
lence of this problem and the solution of a compact matrix equation is pointed out. 

(2) The concise logical formulation of problems involving the rate of approach 
to equilibrium. 


M. R. HEsTENEs in his paper, “Iterative Computational Methods,” is con- 
cerned mainly with matrix problems and with problems that arise from varia- 
tional principles. The main concern is with gradient methods and a discussion of 
the conjugate gradient method of solving linear systems of equations is included. 

J. Topp’s paper entitled ‘‘Motivation for Working in Numerical Analysis,” 
contains a brief discussion of various separate topics in numerical analysis. These 
include evaluation of polynomials, increasing speed of convergence of sequences, 
modified differences, characteristic roots of finite matrices, quadrature, integral 
equations, and game theory. The author stresses the need for an experimental 
and empirical approach to many problems in numerical analysis. 

In the paper, ““Some Numerical Computations in Ordnance Problems,” A. A. 
BENNETT lists and discusses briefly 40 problems considered by the Computation 
Laboratory of the Ballistic Research Laboratories, Aberdeen Proving Ground. 

C. TRUESDELL in the paper, ‘““‘The Simplest Rate Theory of Pure Elasticity,” 
postulates a set of ten non-linear equations for the following ten quantities: The 
density, the three components of the velocity vector field, and the six components 
of the stress tensor. These equations embody the law that the rate of stress is a 
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function of the rate of deformation. Simple special solutions of the equations 
are discussed. 

J. J. STOKER reviews various definitions of stability in his paper, ‘On the 
Stability of Mechanical Systems,” and discusses briefly the possibility of convert- 
ing inherently unstable situations into stable ones by applying periodic external 
forces. 

The longest paper in the collection is the one by F. J. Bureau entitled ‘‘Di- 
vergent Integrals and Partial Differential Equations.’’ In this paper the author 
reviews the definition of the finite part and the logarithmic part of divergent 
integrals and uses these concepts to obtain the solutions to the Cauchy problem 
for the wave equation, the equation for damped wave and the Euler-Poisson- 
Darboux equation. Linear partial differential equations with constant coefficients 
and systems of such equations are also discussed. HADAMARD’s method of singu- 
larities for obtaining solutions to such equations is reviewed and new applications 
of this method are described. The usefulness of the notions of principle and 
logarithmic part of divergent integrals in the theory of elliptic equations is 
illustrated by a discussion of a solution of RADON’s problem. Other related prob- 
lems are also treated. 

A Hh) FT: 


89[1].—E. W. Diyxstra & A. VAN WIJNGAARDEN, Table of Everett's Interpolation 
Coefficients, Computation Department of the Mathematical Centre, Amster- 
dam, Report R294, printed by Excelsior’s Photo-Offset, The Hague, 1955. 
Introduction of 1 page plus 200 pages of tables (pages unnumbered), 24.7 cm. 
The EVERETT interpolation formula for a function f(x) that is tabulated at 

equal intervals of the argument x, = xo + kh is expressible as 


fo = (1 — dD) fo + dfi + Ec? (h)P fo + Ei? (p)H fi + Eo*(p)6*fo 
+ E;*(p)8*f, + Eo®(p)6*fo + Ei®(p)s*fi + ---, 


where fp, = f(xo + ph), fr = f(xx), &"f, are the (2m) central differences of fi, 
and the interpolation coefficients E,2"(p), k = 0, 1, are given by 


E2"(p) = (-y(? » fay '). 


These present tables provide the six functions £,2"(p); k = 0,1; ” = 1, 2, 3; 
for p = 0(0.0001)1; to 7D. Use could have been made of the symmetry relation 
E"(p) = E2"(1 — p) to reduce the size of the book by 50%, either by giving 
E,?"(p) or E;*"(p) alone for ». = 0(0.0001)1 or by giving both functions only as 
far as p = 0.5000. But this present arrangement saves the reader the extra effort 
of looking up the coefficients in two different parts of the book, and it also reduces 
the chances of error. All six functions E,2"(p) are tabulated alongside of each 
other. There are 51 arguments ~ on each page, in groups of five except for the 
last argument on the page which is repeated at the top of the next page. The fine 
interval of 0.0001 in » will, in many cases, obviate the need for interpolating in 
this table. 

These present tables might be compared with two [1] earlier works each 
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containing extensive tabulations of Everett coefficients. A. J. THoMPsON’s excel- 
lent book includes, among other shorter basic tables of £,2"(p), a main table of 
E,?"(p) as far as 2m = 8, for p at intervals of 0.001, with central differences to 
aid in the interpolation for E,2*(p) when p occurs to more than 3D. The NATIONAL 
BUREAU OF STANDARDS volume includes E,*(p) and E;*(p) at intervals of 0.0001 
and Eo*(p), E:*(p), Eo®(p), E:°(p) at intervals of 0.001, all to 10D. In the A-nota- 
tion of this table, EG"(p) = A2%4?(1 — p), Ei"(p) = A2%4?(p). However, the 
user of the Lagrangian table would have the extra work of looking up the desired 
values of £,?"(p) in different parts of the book. Thus in comparison with earlier 
tables this present work has the advantages of finer intervals and more convenient 
arrangement, although it does not have as many decimal places. 

The table was produced on the electronic computer ARRA of the Mathe- 
matisch Centrum in Amsterdam. The program included complete arithmetic 
checks which also covered the signals transmitted to the typewriter. The final 
sheets were reproduced by photo-offset. 

The tables were computed for the ‘‘Empresa nacional bazan de construcciones 
navales militares’ of Madrid.” 

Misprints in the Introduction: line 3, for 2k, read 2n; line 14, for Eo™*(p) 
= E,*(1 — p), read Eo"(p) = E,*(1 — p). 


HERBERT E. SALZER 
Diamond Ordnance Fuze Laboratories 


Washington 25, D. C. 


1. a) A. J. THompson, Tables of the Coefficients of Everett's Central-Difference Interpolation 
Formula. Tracts for Computers, No. V, 2nd ed., Cambridge Univ. Press, 1943. 
b) NationaL Bureau oF STANDARDS, Tables of Lagrangian Interpolation Coefficients, 
Columbia Univ. Press, New York, 1944. 


90[I].—OweEN R. Mock, Fifty Everett integration coefficients. Three typewritten 
leaves (on vellum), deposited in UMT Fine. 


The first fifty-one coefficients to use in EVERETT’s integration formula 


Z1 bed 

[P v@ae = 4D Mes Tyo) + 91, 
=z0 »=0 

listed to 50D and believed by the author to be accurate to 48D. There is no intro- 

duction or description ; the formula is discussed by MILNE [1 ]. The list is believed 

to have been computed on a Model 701 International Business Machine. 


GS. a: 
1. W. E. Mitne, Numerical Calculus, Princeton, 1949, p. 198. 


91[1].—H. E. Sauzer, “‘A simple method for summing certain slowly convergent 
series,” J. Math. Physics, 33, 1955, p. 356-359. 


This gives a table of (the integers) 





B™_;, D®, where A™_; = oye o=(" i ') = BS)_./Dz 
n,"m—t) n ’ n,n—t (m Pe 1)7 i nn, n—t n 


and where D is the least common denominator of the A@_,, ¢=0,1,-:>, 
m — 1. The following values are covered: n = 5,m = 4;n = 10,m = 4,m =7; 
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nm = 15, m = 4, 7, 11; ” = 20, m = 4, 7, 11. The motivation for this tabulation 
is the following: Consider the partial sums s, of a convergent series as a function 
s(x) of x = r. The problem of the determination of s = lim s, is regarded as 
that of interpolation (or rather extrapolation), for s(0) given s(p-') for some 
integral values of p, in particular p = n,n — 1, ---,m — m +1. This is solved 
by Lagrangian interpolation and 


m—1 
$= y nt... 
i=0 


Various standard examples are discussed. The method appears effective when the 
sequence {s,} is monotone. 


io" 


02[1].—H. E. Sauzer, “Equally weighted quadrature formulas over semi-infinite 
intervals,” J. Math. Physics, 34, 1955, p. 54-63. 


In general it is not possible to find real x; = x; such that the quadrature 
formulae 


(1) fC fear = 1/n E feed 
or 
+0 Pal n 
Q) fe pends = = E se 
‘7 M imi 


are exact when f(x) is a polynomial of degree ». The polynomials having the x,“ 
for zeros are obtained formally for » = 2(1)10 in both cases [cf. L. M. MILng- 
THomson, Calculus of Finite Differences, London, 1933, p. 177-180 ]. Among these 
only that for m = 2 in (1) and those for n = 2, 3 in (2) have all roots real. The 
roots are, however, given in all cases to 7D. There follows a detailed discussion 
of a few quadrature formulae of the above type which require exactitude only for 
polynomials of degree m < n, e.g., a seven-point quadrature which is exact for 
cubics in case (1) and a seven-point quadrature which is exact for quintics in 
case (2). 
| A 


93[1].—HerBertT E. Sauzer, Formulas for calculating Fourier coefficients, 5 type- 
written pages (photocopy) placed in the UMT FILe. 

These present formulas enable one to calculate the first twenty-four Fourier 
cosine and sine coefficients of a function f(x) from its values (either given, or 
obtained through calculation or interpolation in a table of f(x) at nine equally 
spaced points ranging from 0 to 27). The coefficients a;, 7 = 1(1)5, and 5, 
4 = 1(1)4, are tabulated for » = 0(1)24, to five decimal accuracy, for use in 
these formulas: 


(1) 2 [5 cos mx dx = a:[f(0) + f(2r)] + as [7(z) +s(7)] 


+a] 4(=) +6(%)] + 0.[s() +7(%)| + een, 
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and 


(2) +f" f(x) sin nx dx = b,[f(0) — f(2x)] + »| 5 (= )- 5(*)| 
f 


+ ols(5) -(F)] +L) - (ZI 


Apart from the error due to the approximate nature of a;, };, and 5("*) ; 


formulas (1) and (2) are exact whenever f(x) is any polynomial of degree <8. 


HERBERT E. SALZER 
Diamond Ordnance Fuze Laboratories 
Washington 25, D. C. 


94[1 ].—HeERBERT E. Saizer, Formulas for inverse osculatory interpolation, 5 type- 
written pages (ozalid) placed in the UMT FILe. 

The author has previously given inverse interpolation formulas [1 ], [2] for 
finding x = x9 + ph from f(x), in terms of f; = f(x;), where the x; = xo + th 
are equally spaced at intervals of h. Those formulas were obtained by the inver- 
sion of LAGRANGE’s interpolation formula. 

The formulas give p in terms of f(xo + ph) = f, f(x:) = fi, and f’(x,) = f’ for 
x; = Xo + that equally spaced intervals h, where the i ranges from — [(m — 1)/2] 
to [n/2], for m = 2(1)7, where n is the number of points required in direct oscula- 
tory interpolation. Although the direct interpolation formula for m = 6 and 
n = 7 is of the 11th and 13th degree accuracy respectively, the inversion formula 
for » which is given does not go beyond the tenth degree terms (and in most 
practical problems it is very rarely that one will go that far; in fact, the first few 
terms will usually suffice). 


HERBERT E. SALZER 
Diamond Ordnance Fuze Laboratories 
Washington, D. C. 


1. H. E. Sauzer, ‘‘A new formula for inverse interpolation,” Amer. Math. Soc., Bull., v. 50, 
er p. 513-516. 
H. E. Sauzer, “Inverse interpolation for eight-, nine-, ten-, and eleven-point direct inter- 
oiladion® J. Math. and Phys., v. 24, 1945, p. 106-108. 


95[K ].—E. S. Pearson & H. O. Hartiey, Biometrika Tables for Statisticians, 
v. 1, Cambridge University Press, 1954, xiv + 238 p. 26.0 cm. Price, $4.50. 
One of the objectives of the originators of the journal Biometrika was to pro- 

vide a place where mathematical tables of importance in statistics could be pub- 

lished. A number of such tables were collated to form the volume I of Tables for 

Statisticians and Biometricians which appeared for the first time in 1914. Four 

subsequent editions were issued and a second volume appeared in 1931. Due to 

the phenomenal growth of new developments in statistical methods since that 
time, it was felt that a modern version was warranted and thus the present 

Biometrika tables. This is the first volume of two and contains fifty-four of the 

more commonly used tables. The second volume will be concerned with the presen- 

tation of more specialized functions. The arrangement of the tables is such as to 
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provide convenience to the user, the grouping being to collect the tabulations of 
individual or related functions. The sources from which the tables were drawn are 
given in an appendix. The older tables, having been in use for a considerable 
time, can be used with the satisfaction of knowing that they are error free. For 
the tables that have been newly compiled one can be certain that painstaking 
efforts, customary from English computers, have been made to guarantee their 
accuracy. 

An outstanding feature of the volume is the comprehensive introduction 
prefacing the tables. Here the underlying mathematical theory is discussed, 
mathematical definitions are provided and a careful exposition is given of the 
method of interpolation to be used in each table, illustrated by example. In 
addition, the usefulness of each table is stressed by the presence of typical appli- 
cations. Teachers of mathematical statistics should find this introduction an 
important supplement to their classroom texts. Because of the excellence of the 
introduction one may consider the volume as a book in mathematical statistics 
supplemented by tables rather than a book of tables preceded by an introduction. 

The tables have been grouped under six main headings. The first group con- 
tains tables of the normal probability function: 


1. The integral P(X) and the ordinate Z(X) in standard units to 7D for 
X = 0(.01)4.5 and to 10D for X = 4.5(.01)6 with first and second central 
differences. This is the fundamental table due to W. F. SHEPPARD. 


2. Values of — log Q(X) = — log[1 — P(X)] to 5D for X = 5(1)50(10) 
100(50)500. 

3. Values of X for extreme values of Q and P to 4D for Q = 0(.0001).02. 

4. Values of X in terms of Q and P to 4D for Q = 0(.001).5. 

5. Values of Z in terms of Q and P to 5D for Q = 0(.001).5. 

6. Table for probit analysis. Tabulates P, Y + Q ® , and a for the 


i . a 2 PQ 
expected probit Y = 1(.01)9. P is given to 3D for Y = .5(.01)6.5, to 4D 
for Y = 6.5(.01)7.5, to 5D for Y = 7.5(.01)8.5, and to 6D for Y = 8.5(.01)9. 
The other functions are given to 4D. 


The second group is concerned with tables derived from the normal proba- 
bility function: 


7. The probability integral of the x? distribution and the cumulative sum of the 
Poisson distribution. Q(x?/v) is given to 5D for x? = .001(.001).01(.01).1(.1) 
2(.2)10(.5)20(1)40(2)134 and » = 1(1) ---. 

8. Values of x? to 5 to 7S (mostly 6S) for Q(x?/rv) = .995, .99, .975, .95, .9, .75, 
5, .25, .1, .05, .025, .01, .005, .001 and » = 1(1)30(10)100. 

9. The probability integral of the ¢-distribution, P(t|v) to 5D for ¢ = 0(.1)8 
and »v = 1(1)24, 30, 40, 60, 120, ©. Upper percentage points of ¢ for 
1 — P(t|v) = .001, .0001, .00001, .000005 and » = 1(1)10 are given to 3 or 
more S. 

10. Chart for determining the power function of the #-test. The probability 6 of 
establishing significance using the two-tail test at the a = 5% or 1% level 
for values ¢ = E(y)/V¥20, is shown on a logarithmic scale. For a = .01, 














11. 


12. 


13. 


14. 


15. 


16. 


17. 


18. 


19. 
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¢ ranges from 2 to 5; for a = .05 from 1.2 to 3.5 for the degrees of freedom 
v = 6(1)10, 12, 15, 20, 30, 60, ©, 6 has the range .1 to .99. 

Test for comparisons involving two variances which must be separately esti- 
mated. The table gives critical values to 2D of the ratio » = (y — 7)/ 
(Aisi? + Aes2”)* for the significance level a = .05 and degrees of freedom, », 
ve = 6, 8, 10, 15, 20, ©; for a = .01, 11, v2 = 10, 12, 15, 20, 30, © and for 
c= AiSi2/ (Ais? + A252?) = 0(.1)1. 

Percentage points of the ¢-distribution. For values of Q = 1 — P(t|v) = .4, 
.25, .1, .05, .025, .01, .005, .0025, .001, .0005, and » = 1(1)30, 40, 60, 120 
and , values of the percentage points are given to 3D. 

Percentage points for the distribution of the correlation coefficient, r, when 
p = 0. The values to 3D are given for Q = .05, .025, .01, .005, .0025, and .0005 
for the degrees of freedom, v = 1(1)20(5)50(10)100. 

Values of z = tanh~'7, with proportional parts are given to 4D for z < 1 
and 3D for z > 1 for r = .000(.002).998. 

Charts giving confidence limits for the population correlation coefficient p, 
given the sample coefficient, r with confidence coefficients .95 and .99. Values 
of r and s are given from —1.(.05)1 for sample sizes = 3(1)8(2)12, 20, 25, 50, 
100, 200 and 400. 

Percentage points of the B-distribution. Values of x are given to 5S for which 


k us (1 — u)>"du 
I.(a, b) = = = .5, .25, .1, 05, .025, .01 and .005, 
f u*(1 — u)>"du 
0 





where a = $v2, 6 = 4», for »; = 1(1)10, 12, 15, 20, 24, 30, 40, 60, 120, 
and vz = 1(1)30, 40, 60, 120, ~. 

Chart for determining the probability levels of the incomplete B-function, 
I,(a, b). Three separate charts are given: Chart A for 1 < 6 < 4, Chart B 
for 4 < b < 15 and Chart C for 15 < 6 < 60. Each chart is further divided 
into upper and lower parts. The lower part contains a family of curves in 
Chart A corresponding to values of a = 2(1)10(2.5)25(5)50(10)100(25)150, 
200. The lower parts of Charts B and C are given at the same intervals except 
that the starting values are a = 5 and 20 respectively. The top part contains 
a family of curves corresponding to selected values of J,(a, b). In charts A 
and B, J.(a, 6) = .005, .01, .02, .05, .1(.1).9, .95, .98, .99, .995, .999, .9995. 
In chart B, in addition, the curve corresponding to the value .9999 is also 
given. In chart C, J,(a, b) = .0005, .001, .005, .01, .025, .05, .1(.1).9, .95, 
.975, .990, .995, .999, .9995. 

Percentage points of the F-distribution (variance ratio). These are given for 
degrees of freedom, »;, = 1(1)10, 12, 15, 20, 24, 30, 40, 60, 120, ©, v», = 1(1)30, 
40, 60, 120, ©, for the upper 25, 10, 5, 2.5, 1, .5 and .1 percentage points to 
2D if F < 100, 1D if F < 1000, to the nearest integer if F < 30,000 and to 
4S for F > 30,000. 

Percentage points of the largest variance ratio s*max/s.?. Values of s*nax/So* for 
the upper 5 and 1% points are given to 2D for the degrees of freedom for 
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$.°, v = 10, 12, 15, 20, 30, 60, ©, for k = 1(1)10. The quantity k here is the 
number of independent variance estimates each based on 1 degree of freedom 
of which s*,,.x is the largest. 


Group three comprises further tables of probability integrals, percentage points, 
etc., of distributions derived from the normal function. 


20. 


21. 


22. 


23. 


24. 


25. 


26. 


27. 


28. 


29. 


30. 


Moment constants of the mean deviation and of the range. For the mean 
deviation this table gives for » = 2(1)20, 30, 60 the expectation to 6D, the 
standard deviation to 4D, variance to 5D and first and second moment ratios 
to 3 or 4D. For the range and the same values of m, values of the expectation 
to 5D, standard deviation to 4D, and variance of the range to 5D, are given 
with the first and second moment ratios to 4 and 3D respectively. The factors 
d,/V, and d,?/V, where d, and V, are the expectation and variance are 
given to 3S. 

Percentage points of the distribution of the mean deviation. The lower and 
upper .1, .5, 1, 2.5, 5 and 10 percentage points of the standardized mean 
deviation as obtained from the probability integral table are given for sample 
size m = 2(1)10 to 3D. 

Percentage points of the distribution of the range. The same lower and upper 
percentage points as in Table 21 of the distribution of the range are given to 
2D for m = 2(1)20. 

Probability integral of the range W, in normal samples of size . The values 
of P(W/n) are tabulated to 4D for m = 2(1)20 and W = 0.00(.05)7.25. 
Percentage points of the extreme standardized deviate from population 
mean (x, — u)/o or (u — x1)/o. The lower and upper .1, .5, 1, 2.5, 5, 10% 
points of X, = (x, — w)/o are given for m = 1(1)30 to 3D. 

Percentage points of the extreme standardized deviate from sample mean 
(x, — £)/o or ( — x:)/c. The lower and upper .1, .5, 1, 2.5, 5, 10% points 
of u = (x, — £)/o are given for m = 1(1)9 to 2D. 

Percentage points of the extreme studentized deviate from sample mean 
(x, — £)/s,or ( — x;)s,. Lower and upper 1 and 5% points of the studentized 
extreme deviate t, = (x, — %)s, where s, is an estimate of o based on » 
degrees of freedom are given form = 3(1)9 andv = 10, 15, 30, © for the lower 
percentage points and = 10(1)20, 24, 30, 40, 60, 120, «© for the upper ones. 
Mean range in normal samples of size m. Values of the expected value, W,, in 
samples from a normal population are given to 5D for m = 2(1)500(10)1000. 
Mean positions of ranked normal deviates (normal order statistics). Values of 
£(i|), the expectation of the i-th ranked normal deviate for m = 2(1)26(2)50, 
4 = 1(1)25 to 3D for m < 20; to 2D for m > 20. 

Percentage points of the studentized range g = (x, — x1)/s,. Lower and upper 
1 and 5% points of g for m = 2(1)20 and » = 10(1)20, 24, 30, 40, 60, 120, 
for lower percentage points and = 1(1)20, 24, 30, 40, 60, 120, © for the 


upper ones. 

Tables for analysis of variance based on range. In part A the scale factor, c, 
to 2D and the equivalent degrees of freedom, », to 1D appropriate to a simple 
classification into K groups of m observations are given for » = 2(1)10 and 














31. 


32. 


33. 


34. 


35. 
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K = 1(i)5, 10, ©, together with constant differences. In part B the like 
values for analysis of double classification with K blocks and n treatments 
for m = 2(1)9 and K = 2(1)10, 20, . 
Percentage points of the ratio s*nax/S*min in a set of K independent mean 
squares. For the upper 5% points, values are given to 3S; for the upper 1% 
points to either 2, 3 or 4S for »y = 2(1)10, 12, 15, 20, 30, 60, © and k = 2(1)12. 
Test for heterogeneity of variance: percentage points of Bartlett's test, M. 
The 5% and 1% significant points of M are given to 2D for k = 3(1)15 and 
c, = 0(.5)5(1)10(2)14. Corresponding to each combination of k and ¢;, two 
percentage points denoted by (a) and (b) are given, representing the maxi- 
mum and minimum values of the true percentage point. 
Test for heterogeneity of variance: table to facilitate interpolation in Table 32. 
This table required only for accurate interpolation in Table 32, gives 
C = c,3/k? and AC = c, — &/k* to 3D for all the marginal entries k and c; of 
Table 32. 
Test for departure from normality. 
A. Percentage points of the distribution of a = (mean deviation) / (standard 
deviation) 
B. Percentage points of the distribution of ¥b; = ms3/mz! 
C. Percentage points of the distribution of bs = m,/m,* 
Part A gives the upper and lower 10, 5 and 1% points for a, to 4D, the expecta- 
tion and standard deviation to 5D for size of sample m = 11(5)51(10)101 
(100)1001. Part B gives the upper and lower 5 and 1% points for ¥b,; to 3D 
and the standard deviation to 4D for size of sample m = 25(5)50(10)100(25) 
200 (50) 1000 (200)2000(500)5000. Part C gives the upper and lower 5 and 
1% points for b, to 2D for m = 200(50)1000(200)2000(500)5000. 
Moments of s/¢ = x/¥Vv and factors for determining confidence limits for o. 
Values of the expectation are given to 6D of the standard deviation to 5D 
and of the standard deviation multiplied by 2V» to 4D. Values of the moment 
ratios 8; = w3*/us* and B, = us/u2? of s/o are given to either 4 or 5D for » 
(degrees of freedom) = 1(1)20(5)100. In addition, the factors ¥»/x_. and 
Vv/x1-2 for confidence limits of ¢ are given for a = .025 and .005 to 3 and 
2D respectively. 


The fourth section gives tables related to various discrete distributions. 


36. 


Test for the significance of the difference between Poisson variables. In 
Part A the values of b are given for which 


6+1 


Eo @Gy sa and 5 MG >e 


for a (the nominal significance level) = .1, .05, .025, .01, .005 andr = a+b 
in the range 1(1)80. Part B for the single-tail test shows for values of 
m = 2.5(2.5)25.0 and the same nominal significance levels as in Part A, the 
probability of rejecting the hypothesis m, = m, when it is true. The prob- 
abilities are given to 3D for a = .1, .05, .025 and to 4D for a = .01, .005. 
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37. 


38. 


39. 


41. 
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Individual terms of certain binomial distributions: 
f(i|n, p) = G)pi(l — p)* 


Values of f(i|, p) are given to 5D for nm = 5,7 = 0(1)5;m = 10,7 = 0(1)10; 
n = 20, 4 = 0(1)20; m = 25, « = 0(1)23; m = 30, « = 0(1)26 and p = .01, 
.02(.02).1(.1).5. 

Significance tests in a 2 X 2 contingency table. The table shows (1) in bold 
type for given a, A and B, the value of b( <a) which is just significant at the 
probability levels of 5, 2.5, 1 and 4% for a single-tail test and (2) in small 
type for the given A, B and r = a + 3, the exact probability to 3D that b 
is equal to or less than the bold type integer. A covers the range from 3(1)15; 
B from 2(1)A, and a from A(—1) until b = 0. 

Individual terms, e~"m‘/z!, of the Poisson distribution. Values are given to 
6D for m = .1(.1)15.0 and z ranging from 0(1) until the contribution becomes 
negligible in the 6th decimal place. 


. Confidence limits for the expectation of a Poisson variable. Tabulates the 


Poisson variable m for the lower and upper .001, .005, .01, .025, .05 confidence 
limits to 3 or 4S for frequencies c = 0(1)50. 

Charts providing confidence limits for p in binomial sampling given a sample 
fraction c/n. Confidence coefficients, .95 and .99. These are charts of the 
confidence belts with the confidence limit P as ordinate and the ratio c/n 
as abscissa for ranging from 8 to 1000. 


Miscellaneous tables such as the PEARSON type curves, rank correlation, and 
orthogonal polynomials are given in the fifth section. 


42. 


43. 


45. 


46. 


47. 


Percentage points of Pearson curves, for given (1, 82, expressed in standardized 
x— pe 





measure. Tabulates the standardized deviate X = to 3S for 6; = 0, 


01, .03, .05, .1, .15, .2(.1)1, and 6, = 1.8(.2)5 for P = .005, .01, .025, .05, 
.95, .975, .99, .995. 

Chart relating the type of Pearson frequency curve to the values of §;, Bo. 
This chart drawn with 8, and {2 as coordinates in a rectangular coordinate 
system associates various types of Pearson curves with points in the 6,62- 
plane and indicates the significantly important regions of application. 


. Distribution of SPEARMAN’s rank correlation coefficients, 7,, in random rank- 


ings. This table gives the probability Q to 3 or 4D of exceeding or equalling a 
given score S,. Results are given for m = 4(1)10, S, = 12(2)308. 
Distribution of KENDALL’s rank correlation coefficient &, in random rankings. 
This table gives the probability Q to 3 or 4D of exceeding or equalling a given 
score for m = 4(1)9, S = 0(2)30; m = 6,7, 10, S = 1(1)35. 

Distribution of the concordance coefficient W, in random rankings. Exact 
values of the probability Q that a given score will be attained or exceeded 
are given for m = 3, m = 3(1)10; m = 4, m = 3(1)6; ” = 5, m = 3. 
Orthogonal polynomials. This table gives values of orthogonal polynomials 
¢:(x) for i = 1(1)6 for sample sizes m = 3(1)52. The observational range 
covered is ¢ = 1(1)m for m = 3(1)12 and the half range ¢ = 1(1)[(m + 1/2) ] 
for n > 12. 
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The final section gives various auxiliary tables which are of fundamental 


importance. 


48. Powers of integers. This table gives the first k powers of the integer n, for 
nm = 1(1)100, k = 2(1)7. 

49. Sums of powers of integers. This table gives the sums of the powers of the 
integers, S(n*) for m = 1(1)100, k = 1(1)7. 

50. Squares of integers. Values of n? are given for m = 0(1)999. 

51. Factorials of integers, their logarithms; square roots; and their reciprocals. 
For m = 1(1)100, 2! and 1/n! are given to 6S, while logy 2!, Vn, 1/¥n, and 
1/n are given to 7D. For m = 101(1)250, m! is given to 6S and logi. m! to 7D. 
For m = 251(1)1000, 7D are given of only logis n! 

52. Miscellaneous functions of » and g = 1 — p over the unit range. For 
p = .00(.01)1.00, (1 + p) and logy {T'(1 + p)} are tabulated to 7D. The 
quantities 1 — p”, pg and p* + g’ are tabulated to4D while ¥1 — p*1, 1 — p* 
and ypq are tabulated to 5D. For p = .910(.002).990(.001)1.000, the quan- 
tities V1 — p°1/V1i — ~ and vVpq are given to 5D. 

53. Natural logarithms, log, x. Log, x is given to 5D together with proportional 
parts for x = 1.00(.05)2.00(.01)9.99. 

54. Useful constants. Mathematical constants, their reciprocals and logarithms 
are given to 10D. Among these constants are multiples and powers of z, e, 
e, logio e, V2, V3, ¥10 and 1 radian. The binomial coefficients (7) are given 
for m = 4(1)12 and i = 0(1)12. Conversion tables from metric to British 
and British to metric are given for measurements of length, area, volume, 
weight and capacity. 








MiILTon ABRAMOWITZ 
National Bureau of Standards 


Washington 25, D. C. 


96[K ].—F. W. Davin & N. L. Jounson, “Statistical treatment of censored data. 
Part I. Fundamental formulae,” Biometrika, v. 41, 1954, p. 228-240. 


This paper is concerned with censored samples in which the smallest k of n 
observations are fully measured, while of the (m — k) largest observations, it is 
known only that x > xz. The authors consider a population described by a con- 
tinuous random variable with p.d.f., f(#), and with distribution function, 
F(X) = S*,f(#dt. Arranged in ascending order of magnitude, observations of 
a (complete) random sample from a population of this type correspond to x,, 
Xo, °**,Xn. Variables X;, Xo, ---,Xz,, defined by the equations F(X,) = r/(n + 1), 
(r = 1, 2, ---, m), are introduced and approximations to the moments of the x,’s 
are found by expanding x, about X, in an inverse TAYLOR series. Central moments 
and cumulants are then obtained from these noncentral moments. With the first 
four moments or cumulants known, functions of the PEARSON system are available 
for approximating distributions of the x,’s and functions of them. 

The first four orders of cumulants of F(x,) are listed in Table 1 which is ap- 
pended to this paper. The first four orders of cumulants of the ordered variables 
%1, X2, ***,%, are listed in Table 2 with an accuracy of O[(m + 2)-*]. Sums of 
powers and products of the first k natural numbers which enter into calculation 
of the cumulants of Tables 1, and 2, are given in Table 3. 
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The application of results obtained in the paper is illustrated by determining 
upper and lower five percentage points of Z,, which is the mean of the first k 
observations (in ascending order of magnitude) of a censored sample from a 
normal population. Calculations are given to 3D for m = 10, k = 2(2)8 and for 
n = 20, k = 4(4)16. 

A. C. CoHEN, Jr. 

University of Georgia 
Athens, Ga. 


97(K ]._F. G. Foster & A. Stuart, “Distribution-free tests in time-series based 
on the breaking of records,’’ Roy. Stat. Soc., Jn., s.B., v. 16, 1954, p. 1-22. 
Two distribution free tests making use of the statistics involving the upper 

records and the lower records for the randomness of a series are introduced in 

this paper. An observation in a time series is called a lower (upper) record if it is 
smaller (greater) than all previous observations in the series. The statistics used 
are s, the sum of the number of upper and lower records in a series and d, the 
difference between the number of upper and lower records. The hypothesis of 
randomness of a series can be tested by the critical region s > sy or |d| > dp. 

The first test provides a test against trend in dispersion while the second tests 

against a trend in location. Comparable round-trip tests (S and D) are also intro- 

duced involving counting the records in the series in both directions. 

The joint distribution of s and d is obtained under the null hypothesis and 
the two variables are shown to be independently distributed as well as being 
asymptotically normal. In the round trip case the corresponding statistics are 
established to be asymptotically normally distributed and the variances have 
been tabulated. 

The powers of the tests against trend in the mean have"been computed using 
empirical sampling techniques on a high speed computer. 

The following tables are given: 


Table 1. P(s < so) to 3 or 4D for m = 3(1)15 and the normal approximation 
for nm = 15. 

Table 2. P(d < do) to 3 or 4D for m = 3(1)6 with normal approximation 
for n = 6. 

Tabie 3. Mean and standard error of s and standard error of d to 3D for 
n = 10(5)100 with approximation for m = 100. 

Table 4. The standard error of D to 2D based on 1000 samples of size n for 
n = 10, 25, (25), 125. 

Table 5. The correlation between d and d’ to 2D for m = 10, 25(25)125. 

Table 6. The power of the d test against normal regression at the 5 percent 
level for trends A = .01(.01).07 and sample sizes m = 25(25)125 based on 1000 
samples in each case. 

Table 7. The power of the D test against normal regression at the 5 percent 
level for trends A = .01(.01).06, .08, and sample size m = 25(25)125, based on 
1000 samples in each case. 

C. F. Kossack 


Purdue University 
Lafayette, Indiana 
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98[K].—H. E. Dantets, “A distribution-free test for regression parameters,” 

Ann. Math. Stat., v. 25, 1954, p. 499-513. 

Consider the usual regression model with y; = ap + Box; + &, where we 
assume only that the x; are distinct and that the errors e; are independent and 
satisfy Pr (e; > 0) = Pr (e; < 0) = 3. The linesa = —x,8 + y,; divide the (8, a) 
plane into polygonal regions, some of them open. Let m denote the minimum 
number of lines which must be crossed to escape from (89, ao) into an open polygon. 
The distribution of m is given to 3D for m = 3(1)30, providing confidence regions 
for (8, ao) and tests for an hypothesis specifying both parameters. The power 
of the m-test is compared with certain competitors, and the complications arising 
when not all x; are distinct are examined. 

J. L. Hopces, Jr. 


University of California 
Berkeley, California 


99[K ].—BENJAMIN EpsTEIn, “Tables for the distribution of the number of 

exceedances,”’ Ann. Math. Stat., v. 25, 1954, p. 762-768. 

Two random samples, each of size m, are taken from the same continuous 
distribution. U,", the number of exceedances, is the number of values in the 
second sample exceeding the rth smallest in the first sample, r fixed. Taking 
advantage of certain symmetries, it is possible to give the complete distribution 
function, Pr (U," < x), in terms of r up to [m/2] andr-—1<x<n-—r-—1. 
Table 1 gives these to not less than 4D and 3S for m = 2(1)15(5)20. [Wayne 
University Technical Report No. 6, by the same author, gives values for 
n = 2(1)20(5)50.] Table 2 indicates, for the case m = 5, how to complete the 
array of probabilities from the portion given in Table 1. 

Leo Katz 


Michigan State University 
East Lansing, Michigan 


100[K ].—C. K. Tsao, “An extension of Massey’s distribution of the maximum 
deviation between two-sample cumulative step functions,” Ann. Math. Stat., 

v. 25, 1955, p. 587-592. 

Let S,(x) and S,,’(x) be the cumulative distributions observed in two random 
samples, of sizes m and m respectively, from the same continuous population. 
Denote by x, and y, the rth smallest in the respective samples, and for 
r < min (m, n) let 

d, = max |Sq(x) — Sn'(x)| 
2<2r 
and 
d= max ' | Sa(x) — Sm’ (x)|. 


x < max (Jr, Xr 
Under the above assumption the distributions of d, and d,’ are independent 
of the distribution of the population. Table I gives P(d, < c/m) to 5D for 


n =m = 3(1)10,15, 20, 30,40; ¢ = 1(1)12; 7 = 2(1)10. Table II gives P(d,’ < c/m) 
to 5D for the same range of arguments. 
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These tables may be used in testing the hypothesis that two populations have 
identical distribution functions. 
FRANK MASSEY 
abrenslty of Oregon 
Eugene, . 


101[K ].—Yosuim1ko Hrraco, Hmenort Mormura, & HisAo WATANABE, 
“Tables for three-sample test,” Inst. Stat. Math., Annals, v. 5, 1954, p. 97-102. 
The authors provide tables for a generalization, to three samples, of a two- 

sample test for randomness originated by WALD and Wo.row!7z [1 ] and studied 

in greater detail by Moon [2]. 

Let 2X1, X2) ***) Xnyy Vis Voy ***» ngs ANA 21, Ze, ***, nz, M1 LS Me < mz, be three 
random samples drawn from the respective populations having continuous 
cdf’s F(x), G(y), H(z). To test whether the m = m, + mz + ms; observations form 
a single random sequence from the same population, i.e., whether F = G = H, 
the values are arranged in a single ascending order and the runs of x, y, z counted, 
regardless of length. Let r;, r2, 73 be the number of runs of x’s, y’s, 2's, respectively 
and put r = 71 +72 +73. For given values of m, m2, ms the critical value of 
r at the a-level of significance is defined as the largest integer, ro, for which 
Prob [r < ro] < a. Table1, for a = .05, is for all m; for which 1 < m, < mz < ms < 10. 
It was obtained from the exact distribution of r discussed by Mood in the above 
paper. Table 2, for a = .05, is for 10 < m, < mz < mz < 30, the n; selected from 
the set 10(2)20(5)30, but excluding ms = m,. Table 3 gives the a = .01 and .05 
critical values for all 2; for which 11 < m, = m2 = ms < 30. Tables 2 and 3 were 
obtained by using the normal approximation given by Mood with a correction 
based on his formulas for the first two moments of 7;. 


P J. LIEBLEIN 
Statistical Engineering Laboratory 
National Bureau of Standards 
Washington, D. C. 


1. A. WaLp & J. WoLrow!7z, ‘“‘On a test whether two samples are from the same population,” 
Ann. Math. Stat., v. 11, 1940, p. 147-162. 
2. A. M. Moon, “The distribution theory of runs,” Ann. Math. Stat., v. 11, 1940, p. 367-392. 


102[K ].—W. J. Drxon, “‘Power under normality of several nonparametric tests,” 
Ann. Math. Stat., v. 25, 1954, p. 610-614. 


Consider a sample of size N, from a normal population with mean py; and 
variance o? and an independent sample of size N2 from a normal population with 
mean ye and variance o”. This paper presents tabulations of the power and power 
efficiency [1] of the rank-sum, maximum deviation, median, and total number 
of runs tests for the null hypothesis yu; = ue with significance level a and alterna- 
tive hypotheses 6 = |: — u2|/o. Table I contains power function values to 3D 
and power efficiencies to 2S for the cases a = 1/10, 1/35, 1/126, Ni = Nz = 3(1)5, 
6 = 0(.25)2(.5)5 for all tests; also for the cases a = 2/126, 4/126, N, = Nz = 5, 
6 = 0(.5)4.5 for the rank sum test; for the cases a = 10/126, N; = Nz, = 5, 
6 = .5(.5)4.5 for the maximum absolute deviation test; and for the cases 
a = 26/126, Ni = Nz = 5, 6 = .5(.5)3 for the median test. Table II contains 


power function values to 3D and power efficiencies to 2S for all tests except 

















_ REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 215 


the total number of runs test for a randomized significance level of a = .025, 
N, = N2=5, 6 = .5(.5)4.5. Table III gives limiting power efficiencies as 
5 — 0 to 3S for the rank sum test for the cases: a = 1/3 and (N;, N2) = (2, 2); 
a = 1/Sand (Ni, N2) = (2, 3);a = 2/15, 4/15 and (Ni, Nz) = (2, 4);a = 1/10, 
1/5 and (Ni, Nz) = (3,3); @ = 1/70 and (Ni, Ns) = (4,4); @ = 1/126, 2/126, 
4/126 and (Ni, Nz) = (5,5);0 < @ < 1 and (Ni, Nz) = (~, @). 

J. E. Wasa 
Military Operations Research Division 
Lockheed Aircraft Corporation 
Burbank, California 


1. W. J. Drxon, “‘Power functions of the sign test and power efficiency for normal alternatives,” 
Ann. Math. Stat., v. 24, 1953, p. 467-473. 


103(K ].—E. S. Pace, “Control charts for the mean of a normal population,” 
Roy. Stat. Soc., Jn. s. B., v. 16, 1954, p. 131-135. 
Consider the well known statistical quality control procedure in which the 
method is to examine a fixed number of items at regular intervals, to measure a 
certain dimension, x, and to plot the sample mean, Z, on a control chart. If the 


: , B , . 
point falls outside the control limits drawn at p + IW (where u is the required 


process mean, and oa is the process standard deviation), the process mean is said 
to be “out of control’’ and some appropriate action is taken. However, it is the 
premise of this paper that unnecessary restriction exists when the width of these 
control limits is fixed, and the author proceeds to present a method by which 
B and N are chosen to obtain a ‘‘best’’ scheme. 

For an estimated process mean, m (perhaps different from yu), and a known 
constant standard deviation, o, it is desired to detect any change in the mean. 
Consider only the case in which we are interested in a single control limit at 
p + Bo/VN and a constant process mean, m. Let P(m) denote the probability 
that a given sample yields a point outside the control limit. Define the ‘average 
run length,”’ L, as the average number of articles inspected between two successive 
occasions when action is taken. The average run length for a constant process 
mean, m, is L = N/P(m). 

The average number of articles produced between two successive occasions 
when action is taken is L/f where f is the fraction of output that is inspected. 
It is easily seen that the average run length is a measure of the amount of scrap 
or defective articles produced when quality is bad, and is also a measure of the 
amount of acceptable articles produced between interferences caused by process 
inspection when the quality is satisfactory. If u is a satisfactory value for the 
mean, m, then for m = yu there is no need for any corrective action to be taken. 
It is desirable that the production should only rarely be suspended. Define Lo as 
the average run length when m = u. Let k (a positive number) be such that a 
shift in the mean of m > u + ko should be detected as soon as possible after it 
has occurred. Define LZ; as the average run length when m = yu + ko. A “‘best”’ 
inspection scheme might be defined as one that minimizes L, for some given large 
value of Lo, or, alternately, as one which maximizes ZL» for some given small 
value of LZ. 
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If only deviations from the mean in one direction are considered, the control 


limit is u + a and Lo , = N/®(B — kvVN), where 


Wand fo = apy?” 
&(x) = (2) f ” exp (—x?/2)dx. 


In accord with the above principles for the one-sided scheme, Tables 1a—c 
give N and B, the latter to 2D, which minimize LZ, for which Ly» = 2,000, 
5,000(5,000)20,000, 40,000, 60,000 is specified for a given k = .2(.1)1.8. Also 
Tables 2a-c give N to 2 or 3S and B to 2D which maximize Ly for which 
L, = 10(5)25(25)100 is specified for a given k = .1(.1)1. 

G. W. McELRATH 
University of Minnesota 
Minneapolis, Minn. 


104[K].—T. E. Sterne, “Some remarks on confidence or fiducial limits,’ 
Biometrika, v. 41, 1954, p. 275-278. 


Define P,,-(7) = (f)a"(1 — x)" and P(a,n,s) = > P,,,-(r) for all r for 


which P,,-(x) < Pn,.(x). From these it is possible to find confidence limits corre- 
sponding to any m, s, and any desired value e, where 7 is the number of trials and 
s the number of successes. It is sufficient to consider the shortest interval of x 
that contains all values of z for which P(z,,s) > ¢ for a given m and s. The 
lower limit is the smallest a satisfying the inequality and the upper limit is the 
largest x. Tables for confidence limits for x to 2 and 3D with confidence coefficients 
0.5 and 0.9 are given for all possible values of s when m = 1(1)10. 


- J. R. VATNSDAL 
State College of Washington 
Pullman, Washington 


105(K ].—H. A. Davin, H. O. Harttey, & E. S. Pearson, ‘The distribution of 
the ratio, in a single normal sample, of range to standard deviation,’”’ Bio- 
metrika, v. 41, 1954, p. 482-493. 

Percentage points of the distribution of w/s, where w is the range and s the 
sample standard deviation of a sample of ” from a normal population, are ob- 
tained partly by the use of an exact distribution applicable in the case of small 
samples, and partly by fitting a curve of the PEARSON system to the first four 
moments for selected values of » and interpolating between these results. The 
final table gives upper a percentage points for m = 3(1)20(10)60(20)100, 150, 
200, 500, 1000 and a = 0.10, 0.05, 0.025, 0.01, 0.005 to 3D for m < 9 and 2D for 
n > 9, and lower percentage points for the same values of a and m > 10. The use 
of this ratio as a test of homogeneity in both large and small samples is considered, 
and it is compared to the ratio of the mean deviation to the sample standard 
deviation and the sample fourth standard moment as a test of normality. 

C. A. BENNETT 


General Electric Company 
Hanford Atomic Products Operation 
Richland , Washington 
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106[K ].—ManrtTIn WEIBULL, “The distribution of t- and F-statistics and of corre- 
lation and regression coefficients in stratified samples from normal populations 
with different means,’ Skandinavisk Akuarietidskrift, v. 36, 1953, supplement 
to haft 1-2. 


This long paper is principally concerned with the application of the standard 
statistics used in exact significant tests in samples from normal to samples from 
stratified universes in which the distributions within each stratum are normal 
with variances constant in all strata but with means different in different strata. 
It begins with a discussion of non-central x?, t- and F-statistics and includes an 
illustrative table of the non-central F distribution in which both the x*’s appearing 
in the numerator and denominator are non-central. 1% and 5% points of this F 
are given to 2D for the degrees of freedom for the numerator, 5; = 1, 2, 4, 8, the 
degrees of freedom for the denominator, b, = 1, 2, 4, 8, 16, 32, ©, but only for 
6b, < be, and the non-centrality parameters 8; and BB, = 0, 3, 1. 


G. J. LIEBERMAN 
Stanford University 


Stanford, Calif. 


107[K ].—A. J. Duncan, ‘‘The use of ranges in comparing variabilities,” Indus- 
trial Quality Control, v. 11, no. 5, 1955, p. 18-22. 


The author tabulates the quantities d.* and » such that E(R/d.*)* = o? and 
the mean and variance of V¥vR/d.*o are those of a x?-distribution with » degrees 
of freedom, where # is the mean of the ranges of k samples of m from a normal 
distribution with variance o*. d,* is given to 2D and » to 1D for k = 1(1)20, 25, 
30, 50 and m = 2(1)7. This is an extension of Table I of Davin [1] where c = d,* 
is given to 2D and » to 1D for k = 1(1)5, 10 and m = 2(1)10, which in turn is an 
extension, with minor corrections, of the original table given by PatTNnarK [2]. 
The methods used for the extension are nec indicated. 


C. A. BENNETT 
Hanford Atomic Products Operation 
General Electric Company 
Richland, Washington 


1. H. A. Davin, “Further applications of range to the analysis of variance,” Biometrika, v. 38, 
1951, p. 393-407. 
2 


. P. B. PatnarK, “‘The use of mean range as an estimator of variance in statistical tests,” 
Biometrika, v. 37, 1950, p. 78-87. 


108[K ].—T. Lewis, ‘99.9 and 0.1% points of the x? distribution,” Biometrika, 
v. 40, 1953, p. 421-426. 


Points of the x? distribution have been computed correct to 8D for the lower 
0.1% points for degrees of freedom m = 1(1)30(10)100, 120, and for the upper 
0.1% points for nm = 40(10)100, 120. 


Using m = 3n,v = 3x7, and P = i™—e-*dt, so that 20 is the 100P% 


I'(m) Je 
point, the author had f(v) = 0 where for n even, 


v iii 
f(v) = tothe +o | = 
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and for odd, 


2 ne v vi yt . 
Se) = {26 fia nat. a’ | aia er || Siete 


Then each value of v was obtained by computing f(v) approximately for an ap- 
proximate value v» of v and then improving by iteration. The first 9) was obtained 
by a formula due to Witson and Hitrerty [1 ], 


Ny oy: ew 7 

m 9m 9m’ 

the x being the normal deviate corresponding to the 100P% value of v. 
Tables of these computations are exhibited. 





T. A. BICKERSTAFF 
University of Mississippi 
University, Miss. 
1. E. B. Witson & M. M. HirFerty, “The distribution of chi-square,” Nat. Acad. Sci., Proc., 
v. 17, 1931, p. 684-688. 


109[K ].—S. H. Aspet-Arty, “Approximate formulae for the percentage points of 

the non-central x? distribution,” Biometrika, v. 41, 1954, p. 538-540. 

The main point of this note is to derive the two approximation formulas of 
the title and compare them numerically with the exact values. The cube root 
transformation of WILSON and HILFERTY turns out to be the appropriate one for 
transforming the non-central x? variate, x”, to an approximately normal variate y. 
The CorNIsH-FISHER expansion [1] is then used to obtain a formula for the 
100a-percentage point of Y (the variate y after standardization) in ascending 
powers of r+, where r = f + \ and f is the number of dégrees of freedom and 
\ is the parameter. This formula is written out to the term in r~!, with the coeffi- 
cients expressed in terms of quantities whose values are given for a = .5(.1).1, 
.05, .01, .005, .001. The numerical comparison with a few exact values indicates 
that even for small degrees of freedom and moderate values of \ the expansion 
agrees very well for both the upper and lower 5-percent points. The probability 
integral approximation for Y is obtained by straightforward application of the 
EDGEWORTH expansion. Numerical comparison with exact values for a few values 
of f and \ show agreement to 3D or better. 


Juttus LIEBLEIN 
Statistical Engineering Laboratory 
National Bureau of Standards 
Washington 25, D. C. 


1. E. A. Cornisa & R. A. FisHer, “Moments and cumulants in the specification of distribu- 
tions,” Extrait de la Revue de I’ Institut International de Statistique, v. 4, 1937, p. 1-14. 


110[K].—A. M. Freupentuat & E. J. GuMBEL, “Minimum life in fatigue,” 
Amer. Stat. Assn., Jn., v. 49, 1954, p. 575-597. 


The authors fit fatigue data by the “survivorship function” 


I(N)s = exp| - (J2teyi]. 
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In this formula, /(N)s is the probability of surviving N cycles under stress S; 
Vz and Nos are location parameters having dimensions of N and the relation 
Vs > Nos > 0; Nos stands for the minimum life measured in cycles to failure 
at the stress level S; 1/as is a dimensionless scale parameter. The problem of 
estimating the three unknown parameters as, No,s, and Vs from a set of fatigue 
data taken at a series of stress levels S is a formidable one. The authors use the 
method of moments to carry out the estimation. This procedure is facilitated by 
the authors’ Table 1. This table gives for 1/as = 0(.01)1(.2)2(1)5 the values to 
4D of the standardized distance from the “characteristic number,” VS to the 
minimum life, of the standardized distance from VS to the mean, and of the third 
central moment in standard units. The method is applied to published data of 
RAVILLy on the fatigue strength under tension of aluminum and nickel wire. 
The question of the sampling errors of estimate obtained by the authors’ pro- 
cedure is left open. 


BENJAMIN EPSTEIN 
Wayne University 
Detroit, Michigan 


111[K].—Max Havperin, SAMUEL W. GREENHOUSE, JEROME CORNFIELD, & 
Juxia ZALokar, ‘“‘Tables of percentage points for the Studentized maximum 
absolute deviate in normal samples,’”’ Am. Stat. Assn., Jn., v. 50, 1955, 
p. 185-195. 


Let xi, 7 = 1 to k, be independent normally distributed variates each with 
mean yu and variance o*. The Studentized maximum absolute deviation is defined 
byd= max ea 4l , where ms*/o? is distributed as x* with m degrees of 

é=1,2,---, k 
freedom and independent of x;. Tables of the upper and lower limit for the upper 
5% and 1% points of d are given to 3S for k = 3(1)10(5)20(10)40, 60 and 
m = 3(1)10(5)20(10)40, 60, 120, o. Examples illustrate the use of the tables. 


L. A. AROIAN 
Hughes Aircraft Co. 
Culver City, Calif. 


112[K].—J. H. CapweE 1, “The probability integral of range for samples from 
a symmetrical unimodal population,” Ann. Math. Stat., v. 25, 1954, p. 803-806. 
An asymptotic expression is given for the probability integral of range for 

samples from a symmetrical unimoda! population. Its accuracy is investigated for 

the case of a normal parent population for sample sizes of 20 to 100. Over this 
range errors are small and by use of a given table of corrections the probability 

integral can be found with a maximum error of 0.0001. The .1, .5, 1, 5, 10, 25, 

50, 75, 90, 95, 99, 99.5, 99.9 percentage points are given to 3S, for sample sizes of 

20(20)100 of the range w from a normal population with unit standard deviation 

and also the mean value of w for the same sample sizes. 


L. A. AROIAN 
Hughes Aircraft Company 
Culver City, California 
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113[K].—S. R. BroapBEnt, “‘The quotient of a rectangular or triangular and a 

general variate,” Biometrika, v. 41, 1954, p. 330-337. 

The paper discusses the distribution of the quotient of the value of a random 
variable from a rectangular or a triangular distribution and the value of a random 
variable from an independent distribution. The denominator distribution is 
restricted to one for which a frequency function exists almost everywhere and for 
which, if the cumulative distribution function is F(x) = z, the inverse function, 
G(z) = x, is defined, non-zero, and finite for almost all z. 

The tables of the paper are for the cases in which the denominator distribution 
is a ‘‘normal” distribution. Table 1 is entitled ‘‘Percentage points of the quotient 
of a rectangular and an independent normal variate. If y is rectangular with range 
(1 — a, 1+ a) and x is normal with mean 1 and variance §?, the table gives 
percentage points of (y/x).’’ The values given are solutions for Q to 3D of [1] 


1 BQ 1-—a 1 i+a 1 
me(—1) + 8c omg) —m(52-2)f = 
hl, Ql. \202...8/ “Le. #8 
k = .01, .05, .95, .99; a = 0(.02).10; 8 = 0(.01).05. Table 2 is entitled ‘‘Per- 
centage points of the quotient of a triangular and an independent normal variate. 
If y is triangular with range (1 — 2a, 1 + 2a), and x is normal with mean 1 and 


variance 6’, the table gives percentage points of (y/x).’’ The values given are 
solutions for Q to 3D of 


mu $) ae (S23) - (3-2) 


1 + 2a 1 
m(!22*—1)] 
as: Se ee 














k = .01, .05, .95, .99; @ = .00(.01).05; 8 = .00(.01).05. 
In all cases tabled, Hho(—1/8) is negligible and could be omitted from the 
above formulas without affecting the solution. 


K. J. ARNOLD 
Michigan State University 
East Lansing, Mich. 


1. For the definition of Hh,(x) see BAAS Math. Tables, v. 1, London, 1941, p. x. 


114[K].—S. B. Cnowpuury, “The most powerful unbiased critical regions and 
the shortest unbiased confidence intervals associated with the distribution of 
the classical D*-statistic,” Sankhyd, v. 14, 1954, p. 71-80. 

This continuation of an earlier paper by P. K. Bose [1] includes two short 
tables for determining the shortest unbiased upper and lower 95% confidence 
intervals for a parameter A? used in testing the equality of means of two p-variate 
normal populations. Let the (known) variance matrices of the two populations 
be A = (ai;), A’ = (a’;;), and their unknown vector means be {a;}, {a}. In 
terms of these quantities, 


A? = pL VBi(a: — ai’) (a; — a;'), 
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where 64 are the elements of the inverse matrix (8;;)—, Bij = (n’ag-+-na’ 4) /(n+n’). 
To test the equality of means a random sample is drawn from each population, 
of, say, respective sizes n, m’, and sample means {a;}, {a,/}. One then computes 
a sample statistic D? (to which the “classical D®”’ statistic introduced by 
MAHALANOBIs is related by D? = D;? — 2/f), corresponding to A’, by means of 
D? = p“L YB (a; — a’) (a; — a;’), and from this, the quantities P = kD,*/2, 
L = (a)4l, where i is the harmonic mean of n, n’. 

Table 1 gives, for p = 1(1)10 and for 14 values of L ranging from L = 2 to 
L = 400, the lower 95% confidence limit \ for the parameter \ related to A? in 
the same way as / is related to D,’, for all alternatives \ > Ao. Here A» represents 
the hypothesis to be tested, that is, A» represents the bound for A which would 
assure that A? is sufficiently small to make the means equal within desirable 
tolerances. The values of \ = \(L) are the solutions, for the given values of L, of 


[ove amar = 1 - a, 


where a = .95, and 
V(L, d) = L?)-@-D 4-H HONT 0 LX) 


(where J is the imaginary Bessel function) is the probability density function of L. 

Table 2 gives, for the same values of p and L, the upper 95% confidence 
limit \, for all alternatives 4 < Xo, in precisely analogous manner to Table 1, 
with the limits of integration 0 to L in the equation for \(Z) being replaced 
by L to ~, 


Jutrus LIEBLEIN 
Statistical Engineering Laboratory 
National Bureau of Standards 
Washington 25, D. C. 


1. P. K. Bosg, “On recursion formulae, tables, and Bessel function populations, associated 
with the distribution of classical D*-statistic,”” Sankhyd, v. 8, 1947, p. 235-248. 


115[K].—D. R. Cox, “‘The mean and coefficient of variation of range in small 
samples from non-normal populations,”’ Biometrika, v. 41, 1954, p. 469-481. 
The author is interested in the mean range and the coefficient of variation 

of the range in small samples of 2, 3, 4, and 5, from non-normal populations. 

Different types of populations covering a wide range of values of 8; and B2, the 

usual measures of skewness and kurtosis, are considered: symmetrical and un- 

symmetrical mixed normal distributions, the normal distribution, the rectangular 
distribution, exponential type distributions, the PEARSON system, and SHONE’s [1 ] 
numerical results for 5 populations of discrete values. Based on these results he 
provides a table for the normalized mean range and the coefficient of variation of 
the range to 3D for sample sizes of 2, 3, 4 and 5, for B, = 1(.2)2(.5)5(1)9; 8B: is 
not a determining factor. Comparisons are made between the distribution function 
of the range and the ratio of two ranges from the exponential e~* and the unit 
normal population. The theory is applied to the point estimation of dispersion 
by use of the range, the control chart for the range, the comparison of several 
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ranges by use of the range ratio test and the use of the range in the ¢-test for non- 
normal as well as normal distributions. Many short tables give numerical illustra- 
tions. The author concludes that if 82 is known, it is possible to make a rough 
correction for non-normality in methods that use ranges of small samples, and 
that on the whole these methods are less affected by non-normality than corre- 
sponding methods using variances of small samples. 


L. A. AROIAN 
Hughes Aircraft Company 


Culver City, Calif. 


1. K. J. SHonE, ‘Relations between the standard deviation and the distribution of range in 
non-normal populations,” Roy. Stat. Soc., Jz., s.B., v. 11, 1949, p. 85-88. 


116[K ].—H. A. Davin, “The distribution of the range in certain non-normal 

populations,” Biometrika, v. 41, 1954, p. 463-468. 

Explicit expressions for the probability integral P(w|m) and the expectation 
E(w|n) of the range w of a sample of m are obtained for several non-normal 
distributions with varying degrees of skewness and/or kurtosis. Values of 
Q(w|n) = 1 — P(w|n) for w = k(n, a)oz, where k(n, a) is the 100a% point of 
w/o for the normal distribution, are given to 3D in Table 1 for a = .01, .05, 
nm = 2(2)12, and various degrees of skewness in the populations considered. Part 
of these values were obtained from the explicit expressions and part by direct 
quadrature using related tables. Table 2 gives the ratio E(w|m)/d,o, to 3D, where 
d, is the expected range, in standard units, of a sample of m from the normal 
distribution, for m = 2(2)12 and those distributions for which E(w|m) was ob- 
tained explicitly. Table 2 indicates that appreciable non-normality has little 
effect on the use of the range to estimate the standard deviation, but Table 1 
further substantiates the difficulty in distinguishing between changes in varia- 
bility and departures from normality. 


C. A. BENNETT 
General Electric Company 


Hanford Atomic Products Operation 
Richland, Michigan 


117[K, L].—H. J. Gawitk, A table of a function related to the error function, 
Armament Research and Development Establishment, BR. Memo. B4/1/55, 
Fort Halstead, Sevenoaks, Kent. 3 typescript pages mimeographed, 10 leaves 
of tables (photo-offset) stapled, 33 cm. Deposited in the UMT Fine. 


The table lists F(x) to 10D for x = 0(.01)10.00, where 


F(x) = e#*? f e”/*dt, 
0 
There is a brief discussion of the computational techniques used and the 
author’s estimate that all digits are accurate except for possible round off in- 
accuracy in the last digit. There is no discussion of use of the table or interpolation. 
The calculation was on the AMOS machine. 


C; BF. 
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118[L ].—OweEn R. Mock, Table of Bessel functions at multiples of x/2. 8 leaves 

(mimeographed both sides), deposited in UMT Fixe. 

A table of J,(x) for m = 0(1)59 and x = kx/2 to 20D for all non-negative 
integral values of k giving significant values at this precision and not exceeding 25. 
There is no introduction or description. The tables are believed to have been com- 
puted on a Model 701 International Business Machine, and the author believes 
that all digits are correct. 

Cc. B..F. 


119[L ].—ALBert D. WHEELON & Joun T. Ropacker, A Table of Integrals 
involving Bessel Functions, The Ramo-Wooldridge Corporation, Los Angeles, 
California, 1954, 91 p., 28 cm. Multilithed from manuscript and typescript. 
Deposited in UMT FILe. 

The authors note that the preparation of this table of 400 integrals was halted 
when they learned of the BATEMAN Manuscript Project which was doing similar 
work. They hope the tables contain the most frequently encountered integrals. 
This set is much less comprehensive than those in the Bateman Manuscript 
Project, ‘“Tables of Integral Transforms,” v. 2, New York, 1954, p. 5-177. There 
are, however, a few formulas in the present work which are not found in ‘“Tables 
of Integral Transforms.” 

The authors’ careful indication of sources of their formulas is valuable. 

Sometimes the statement of the formula is careless. In particular, the reviewer 
noted the following: 


p. 40, 11.306. The important condition b < a is not stressed. 

p. 41, 11.311. The interpretation of the formula for a > 1 might be noted. 

p. 46, 11.201. Many similar integrals published by W. N. BatLey, “Some 
infinite integrals involving Bessel functions,” London Math. 
Soc., Proc., v. 40, 1936, p. 37-48, are not listed. However, this 
particular integral seems new. 

p. 15, 1.207. The integral on the left is divergent. The authors probably 
intend to have the principal value used. 


The indexing scheme is clearly described, and integrals are easily found 
in the list. 


P. HENRICI 
National Bureau of Standards 
Washington, D. C. 


120[M, S].—A. J. W. Dutjvestijn & J. BerGuuts, The computation and the 
expansion of some triple integrals originating from the theory of cosmic rays. 

R 261, Math. Centrum, Amsterdam. Seventeen leaves mimeographed from 

typescript deposited in the UMT FILE. 

This work contains tables of K and J(r) for a = 0.005, 0.01, 0.05, 0.1, 0.5, 1 
and J*(B) for a = 0.1 and 1. Two cases are considered: (1) a = 0.0245/r', 
B = 0.0245/r:2, and y = 0.0245/r.; (2) Ina = In 0.454 — Inr + In (1 + 4r) 
— 47°83, In B = In 0.454 — Inr,; + In (1 + 4r;) — 47°", In y = In 0.454 — Inn; 
+ In (1 + 4r2) — 4r,?/. 
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Here 71,2 = (a? + r? + 2ar cos@)"?, and 
I(r, B) =r es (1 — ee) (1 — €-8*)(1 — 870 
J*(B) = B-8P £ I(r, B)dr 
I(r) = f ” B-SA](r, B)dB 
_ 


K= J(r)dr = f J*(B)dB. 
r=0 B=0 
Behavior of J*(B) for small and large positive values of B is discussed, and 
the method of computation is described. 
The integrals are of interest in studies of cosmic rays. 
eth Ay a 


TABLE ERRATA 


In this issue references have been made to errata in Review 80, Review 81, 
and Review 89. 

In MTAC, 9, July 1955, reference was made to errata in Review 55 and 
Review 59. 


245.—E. JAHNKE & F. Empe, Tables of Functions. Dover edition, 1945. 


On p. 211, section 9, chapter VIII, the left side of the last formula should read 
— 2,(z) + N,(z) --- instead of 2,(z) + N,(z) ---. 


Henry E. FEttI1s 
931 Five Oaks Avenue 
Dayton 6, Ohio 


246.—ROBERT E. GREENWOOD, ‘“‘Coupon Collector’s Test for Random Digits,” 
MTAC, 9, 1955, p. 3. 
Entry for m = 13 should read .00800 8315 2 instead of .0080 9315 2. 


ROBERT E. GREENWOOD 
The University of Texas 
Austin, Texas 


NOTES 
R. C. Archibald 


Professor RAYMOND C. ARCHIBALD, who founded Mathematical Tables and 
Other Aids to Computation and who was always one of its principal contributors and 
the most valued adviser to the Editorial Committee, died on July 26, 1955 at 
Sackville, New Brunswick, Canada. 

Professor Archibald was the Chairman of the Editorial Committee from the 
first issue of the journal in 1943 through 1949. 

A short history of some of Professor Archibald’s professional and personal 
accomplishments will appear in MTAC shortly. 
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There can be no adequate statement of the loss felt by the Editorial Com- 
mittee and other people interested in mathematical tables, computation, and 
related topics on learning of Professor Archibald’s death. 

Ge Bt. 


Postdoctoral Research Associates of the 
National Bureau of Standards 


The National Academy of Sciences, the National Research Council, and the 
National Bureau of Standards announced a plan for postdoctoral research asso- 
ciates at the National Bureau of Standards in October 1954. The program is 
intended to provide advanced training, through Bureau sponsored research and 
through close association with outstanding members of the staff of the Bureau, 
to a few carefully selected young scientists. 

This program is directed by Dr. WALLACE R. Brope, Mr. JosEpH HILsENRATH, 
and Dr. Davip E. Mann of the National Bureau of Standards. 

Three appointments for the current year are in numerical analysis or closely 
related fields. 

Dr. RICHARD J. PROSEN plans to apply the computer SEAC to an investigation 
of crystal structure, with particular emphasis on problems of atomic arrangement. 
He hopes to answer the question, “How much and how well can an electronic 
computer be used in such problems?” (Mr. Howarp R. McMurpre, NBS 
advisor.) 

Dr. Morris Krauss plans a theoretical quantum mechanical study of the 
relationship between electron energy levels and molecular structure. Such calcula- 
tions may help explain certain basic problems that are important in other fields 
such as mass spectrometry. SEAC will be the primary tool used in the study. Dr. 
Krauss will attempt to answer two questions: ‘“‘How do you go about solving 
such a highly complex problem?’ and “Of what value are known methods?” 


National Simulation Conference 


The Dallas-Fort Worth Chapter of the Institute of Radio Engineers Profes- 
sional Group on Electronic Computers (PGEC) will sponsor a National Simula- 
tion Conference in Dallas, Texas, on 19-21 January 1956. 

The Conference will be devoted to simulation and associated computing 
techniques, and will include topics in (a) general simulation (mathematical, 
physical, logistic, etc.) ; (b) advances in computer design, techniques, and applica- 
tions; and (c) methods of determining and improving the accuracy of analog 
solutions. 

While the I.R.E. does not sponsor meetings requiring military security clear- 
ance, an attempt will be made to arrange such a session on an “‘associated”’ basis 
if interest warrants. 

It is planned to publish “proceedings” of the Conference. 

Further information regarding the Conference can be obtained from 


Mr. J. R. Forester 
2104 Huntington 
Arlington, Texas. 
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The Second Regional Graduate Summer Session in Statistics 
held at the University of Florida 
in June and July 1955 


At the request of the Southern Regional Education Board's Advisory Com- 
mission on Statistics, the University of Florida, Virginia Polytechnic Institute 
and North Carolina State College agreed to initiate a continuing program of 
graduate summer sessions in statistics to be held at each of the three institutions 
in rotation, beginning in the summer of 1954. At the first session held at Virginia 
Polytechnic Institute eighty-nine students were enrolled, and at the second session 
ninety-six students were enrolled. The third session will be held at North Carolina 
State College in Raleigh in the summer of 1956. 

It is the purpose of this program to serve: (1) research and professional workers 
who want intensive instruction in basic statistical concepts and modern statistical 
methodology; (2) teachers of elementary statistical courses who want formal 
training in modern statistics; (3) graduate students in other fields who desire 
supporting work in statistics; (4) prospective candidates for graduate degrees in 
statistics, and (5) professional statisticians who wish to keep informed about 
advanced specialized theory and methods. 

Each of the summer sessions will last six weeks and each course will carry 
approximately three semester hours of graduate credit. 


Sponsorship by the Office of Naval Research, the National Bureau 
of Standards, the National Science Foundation, 
and various other organizations 


Much of the editorial and refereeing work for this and for earlier volumes of 
Mathematical Tables and Other Aids to Computation has been done under the 
sponsorship of the Office of Naval Research and the National Bureau of Standards 
as part of their program for the development of numerical analysis and research 
connected with numerical analysis. Considerable help has been given to the 
Editorial Committee by various committees sponsored by the National Science 
Foundation and several learned societies in connection with mathematical and 
statistical tables. Finally, the members of the Editorial Committee who are 
members of University faculties and the referees of papers submitted to the 
journal have carried out their work for Mathematical Tables and Other Aids to 
Computation as part of their University research assignment. The contribu- 
tions of all these organizations to the publication of this journal are gratefully 


acknowledged. 


Index 


In this issue the Editorial Committee has decided to try to set a format for a 
more useful index of reviews and other notices of literature available. In particular, 
we try here to furnish a classified list of tables which will help readers who wish 
to use Mathematical Tables and Other Aids to Computation as a supplement to 
existing indices of tables, or who are looking for some particular table. 
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It is not feasible to repeat enough of the review to make the index compete 
with such works as ‘‘An Index to Mathematical Tables’’ by Fletcher, Miller, and 
Rosenhead ; however, in order to make the index as convenient as is feasible all 
entries are given by full title as published in the review indexed. In the review the 
reader is likely to find a description of the table including: 


(1) the precision (a number followed by D indicates the number of digits 
to the right of the decimal point listed, a number followed by S indicates the 
number of significant decimal digits listed) ; 

(2) the range of the arguments and their increments (in a standard nota- 
tion the first number listed is the smallest value of the argument, the standard 
increment following this smallest value is listed in parentheses, then the largest 
value attained by steps of this increment is listed, then the next increment is 
listed in parentheses, and so on) ; 


(3) other pertinent information including uses, related tables and works, 
and so on. 


The index of reviews will not be limited to tables. It is divided into several 
principal sections in the following order (in this volume): 


General Works on Numerical Analysis 
Number Theory 

Algebra 

Analysis 

Statistics 

Physical Sciences 

Automatic Computation 


The entries under Analysis are the most voluminous. These include all the 
mathematical tables which would normally be listed in ‘‘An Index of Mathe- 
matical Tables’’ by A. Fletcher, J. C. P. Miller, and L. Rosenhead, Scientific 
Computing Service, London, 1946, except those in number theory or algebra. In 
addition, some tables whose entries are restricted by physically measured quan- 
tities but which otherwise are analytical in nature may be included. 

The classification of the tables indexed under Analysis has followed that of 
the work cited above through the first digit to the right of the decimal point. 
No further subdivision has been attempted. In a few cases new subdivisions have 
been added where no appropriate subdivision was found to exist in the earlier 
work. The sections are in the same order as the ones listed in the cited work, but 
the numerical headings have been omitted. The verbal headings of the sections in 
this work have been retained more or less intact; however, they have been ex- 
panded in some cases so that they seem to be self-sufficient without the necessity 
of including the whole indexing scheme. The Editorial Committee has corre- 
sponded with Messrs. Fletcher, Miller, and Rosenhead about this index and we 
are pleased to acknowledge the cooperation of these authors in our indexing work. 

The only major deviations from the usage in the Index of Fletcher, Miller, 
and Rosenhead have resulted from the recent popularity of numbers of non- 
decimal radix, particularly in connection with automatic computers, and the 
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addition to the index of approximations, such as the ‘“‘Approximations for Digital 
Computers” by C. Hastings, Jr., J. T. Hayward, and J. P. Wong, Jr., Princeton, 
1955. References to numbers expressed in common integral bases (such as base 
two) or expressed symbolically without reference to base (such as some of the 
approximations cited above) have been made under headings which might have 
been restricted to decimal headings in the Imdex copied. However, logarithms to 
bases which are non-decimal have been assigned a special section, not used by 
Fletcher, Miller, and Rosenhead, and immediately following the section dealing 
with the corresponding decimal logarithm. 

The classification of the tables relating to Number Theory has been carried 
out in the same way with reference to the “Guide to Tables in the Theory of 
Numbers” by D. H. Lehmer, The National Research Council, Washington, D. C., 
1941. The full detailed classification employed by Lehmer has been used here. 
Again, only the verbal headings have been repeated, with some augmentation 
when this seemed desirable to make them stand alone with meaning, and the 
numbering symbols have been omitted. The order of the listings is the same as 
that of Lehmer. The Editorial Committee acknowledges the cooperation of 
Professor Lehmer in granting permission to use his classification. 

The classification of tables in Statistics has been done largely by Dr. Bernard 
Sherman using a complete outline kindly furnished to the Editorial Committee 
by Professor H. O. Hartley, who is using it in his work on indexing Statistical 
Tables; we are pleased to acknowledge his help in making his index headings 
available to us. Tables which might be classified in the Algebra or Analysis sec- 
tions (some tables which might relate to block designs of experiments and some 
tables of functions which—like the hypergeometric function—are analytic in 
origin but have considerable applications in statistics) may not always be in- 
cluded in the Statistics index; however, no attempt has been made to avoid 
multiple listings, and such omissions are oversights. Tables of quality control 
have been listed under the miscellaneous heading of the normal distribution. 

Adequate listings for the other classifications seem not to be obtainable. At 
the moment, the number of entries is sufficiently sparse to lead to a belief that 
the reader will have little difficulty scanning the headings and finding those which 
seem likely to hold material of interest to him in his particular search. 

It should be noted that the headings used are not the only ones which will 
be used in later indices. In particular, it seems likely that a new category of Social 
Sciences will be introduced in the next volume. 

With regard to many of the more elementary functions, it should be noted 
that the policy of MTAC does not include review of tables of the most common 
functions. Hence, those included in this index are only those found in more 
comprehensive books of tables reviewed, and the sections of the index relating to 
the most common functions should not be considered as even approximately 
indicating the range of publications including values of these functions during the - 
period reviewed. 

The Editorial Committee hopes to be able to assemble a complete index of 
reviews for the first ten volumes of MTAC shortly after the tenth volume has 
been completed. It will be appreciated if authors of works reviewed and others 
will suggest proper categories for these works without further solicitation. Com- 
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ments on the utility of this index and corrections to the present listing will be 
appreciated if they are sent to the Chairman of the Editorial Committee. 
Except for the help cited and considerable advice from Drs. T. H. Southard 
and M. Melkanoff, the index has been arranged completely by the Chairman of 
the Editorial Committee, for the first attempt seemed necessarily a job to be 
undertaken in one place. Mistakes and errors of judgment are wholly his. 
me Be 
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RoBERT E. GREENWOOD, “Coupon Collector’s Test for Random Digits,” 
MTAC, v. 9, p. 1, the displayed relation for p, should read 


=. {9 — 
Pn = oa = (—1)/ (’) (9 Be D) ’ 
instead of 


i (4 = 
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E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 
Algebraic Constants. Roots of Algebraic Equations 
ANDRES ZAVROTSKY Tablas para la Formacion de las Ecuaciones Cubicas 86 198-199 
Conversion Tables 
FRIEDRICH Léscu Siebenstellige Tafeln der elementaren transzendenten 78 195-196 
Funktionen 
E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 


H. O. HaRTLEY 


Fundamental Tables of Logarithms to Non-decimal Integral Bases 


Cecit HASTINGS, JR., Approximations for Digital Computers 56 121 
JEANNE T. HAYWARD, 
James P. Wonca, JR. 


Natural Trigonometricai Functions with Argument in Radian 


NBS Applied Tables of Circular and Hyperbolic Sines and Cosines 33 71 

Mathematics for Radian Arguments 

Series No. 36 

W. FLicce Four-place Tables of Transcendental Functions 57 123 

NBS Applied Tables of Sines and Cosines for Radian Arguments 58 124 

Mathematics 

Series No. 43 

Frrepricu Léscu Siebenstellige Tafeln der elementaren transzendenten 78 195-195 
Funktionen 


Natural Trigonometrical Functions with Argument 
in Degrees and Decimals 


W. FLtcce Four-place Tables of Transcendental Functions 57 123 


Natural Trigonometrical Functions with Centesimal Argument 


Cecit HASTINGS, JR., Approximations for Digital Computers 56 121 
JEANNE T. HAYWARD, 
James P. Wonca, Jr. 


Inverse Circular Functions 
W. L. STEVENS “Tables of the angular transformation” 6 28 
Cecit HAsTINGs, JR., Approximations for Digital Computers 56 121 


JEANNE T. HAYWARD, 
James P. Wona, Jr. 


FRIEDRICH Léscu Siebenstellige Tafeln der elementaren transzendenten 78 195-196 
Funktionen 











Author /s 


W. FLUGcEe 


FRIEDRICH Léscu 


E. S. PEARSON 
H. O. Hartley 


TABLES FROM ANALYSIS 


Title 
Natural Values of e* 
Four-place Tables of Transcendental Functions 


Siebenstellige Tafeln der elementaren transzendenten 
Funktionen 


Biometrika Tables for Statisticians 


Quantities (Natural and Logarithmic) connected with e 


Ceci Hastincs, JR., 
JEANNE T. HAYWARD, 
James P. Won, Jr. 


W. FLUcce 


NBS 
Appl. Math. 
Ser. No. 36 


W. FLUGGE 


FrrepricH Léscu 


Cecit HAsTINGs, JR., 
JEANNE T. HAYWARD, 
James P. Wonca, Jr. 


W. FLUcce 
FRIEDERICH LéscH 


E. S. PEARSON 
H. O. HartLEy 


Approximations for Digital Computers 


Four-place Tables of Transcendental Functions 


Natural Values of Hyperbolic Functions 


Tables of Circular and Hyperbolic Sines for Radian 
Arguments 


Four-place Tables of Transcendental Functions 
Siebenstellige Tafeln der elementaren transzendenten 
Funktionen 

Natural Logarithms of Numbers 
Approximations for Digital Computers 


Four-place Tables of Transcendental Functions 


Siebenstellige Tafeln der elementaren transzendenten 
Funktionen 


Biometrika Tables for Statisticians 


Natural Logarithms of Reciprocals, Square Roots and other Functions 


Ceci HASTINGS, JR., 
JEANNE T. HAYWARD, 
James P. Wona, Jr. 


FRIEDERICH Liscu 
E. S. PEARSON 


H. O. HARTLEY 


AECD-3497 
S. FRANKEL, 
E. NELSON 


W. Fitcce 


Approximations for Digital Computers 


Inverse Hyperbolic Functions 


Siebenstellige Tafeln der elementaren transzendenten 
Funktionen 


Biometrika Tables for Statisticians 


Miscellaneous Functions of Circular and Hyperbolic 


Functions and Inverse Functions 


“Methods of treatment of displacement integral 
equations” 


Exponential Integrals Ei (x) 
Four-place Tables of Transcendental Functions 





Review 

Number Page 
57 123 
78 195-196 
95 205-211 
56 121 
57 123 
33 71 
57 123 
78 195-196 
56 121 
57 123 
78 195-196 
95 205-211 
56 121 
78 195-196 
95 205-211 
35 72 
57 123 











TABLES FROM ANALYSIS 239 


Review 
Author /s Title Number Page 
Exponential Integrals. —Ei(—x) 
Cecit HAsTINGs, JR., Approximations for Digital Computers 56 121 
JEANNE T. HAYWARD 
James P. Won, Jr. 
Sine Integral 
AKADEMIIA Nauk SSSR Tablitsy Integralnogo sinusa i kosinusa [Tables of 39 74 
the sine and cosine integral] 
NBS Applied Table of Sine and Cosine Integrals for Arguments 23 37 
Mathematics from 10 to 100 
Series No. 32 
W. FLicce Four-place Tables of Transcendental Functions 57 123 
Cosine Integral 
NBS Applied Table of Sine and Cosine Integrals for Arguments 23 37 
Mathematics from 10 to 100 
Series No. 32 
AKapemiiA Nauk SSSR Tablitsy Integralnogo sinusa i kosinusa [Tables of 39 74 
the sine and cosine integral} 
W. FLUGGE Four-place Tables of Transcendental Functions 57 123 


The Gamma Function 


Cecit HastTincs, JR., Approximations for Digital Computers 56 121 
JEANNE T. HAYWARD, 
James P. Wonca, Jr. 


W. FLicce Four-place Tables of Transcendental Functions 57 123 
E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 

Logio (1 + x) 


E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 

Gamma Function for Complex Arguments 
NBS Table of the Gamma Function for Complex Arguments 38 74 
Applied Math. 
Series No. 34 

Incomplete Beta Function 

E. S. PEARSON Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 

The Error Integral. Tables of the Ordinate 


(the Integrand) 

Cecit HASTINGS, JR., Approximations for Digital Computers 56 121 
JEANNE T. HAYWARD, 
James P. Won, Jr. 
E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 

Tables of the Error Integral 
Cecit HastIncs, JR., Approximations for Digital Computers 56 121 


JEANNE T. HAYWARD, 
James P. Won, Jr. 











Author /s 


W. FLtcce 


E. S. PEARSON, 
H. O. HartTLEy 


E. S. PEARSON, 
H. O. HARTLEY 


H. J. GAWLIx 


Ceci Hastincs, JR., 


TABLES FROM ANALYSIS 


Title 
Four-place Tables of Transcendental Functions 
Biometrika Tables for Statisticians 


The Error Integral. Ratio of Integral to Ordinate 
Biometrika Tables for Statisticians 
A Table of a Function related to the Error Function 


The Error Integral. Inverse Tables 
Approximations for Digital Computers 


JEANNE T. HAYWARD, 
James P. Wonca, Jr. 


E. S. PEARSON, 
H. O. HARTLEY 


E. S. PEARSON, 
H. O. HARTLEY 


K. G. TapMAN 


CHESTER SNOW 


L. Fox 


W. FLUcce 
OweEn R. Mock 


L. Fox 


W. FLUcce 


L. Fox 


W. FLtcce 


Biometrika Tables for Statisticians 


Moments. Chi-squared Tables 
Biometrika Tables for Statisticians 


The Error Integral with Complex Argument 
“‘Hohenkarte des Fehlerintegrals” 


Legendre Polynomials and their Applications 

Tables of Flow Functions for Bodies of Revolution in 
Circular Tunnels and Jets 

Derivatives of Legendre Polynomials and their Applications 
Tables of Flow Functions for Bodies of Revolytion in 
Circular Tunnels and Jets 

Legendre and Associated Functions of Fractional 
or Complex Order 
“Formulas for computing capacitance and in- 
ductance”’ 
Bessel Functions of the First Kind: Integral Order 


A Short Table for Bessel Functions of Integer Order 
and Large Arguments 


Four-place Tables of Transcendental Functions 
Table of Bessel Functions at Multiples of x/2 


Bessel Functions of the Second Kind: Integral Order 


A Short Table for Bessel Functions of Integer Order 
and Large Arguments 


Four-place Tables of Transcendental Functions 


Modified Bessel Functions of the First Kind: Integral Order 


A Short Table for Bessel Functions of Integer Order 
and Large Arguments 


Four-place Tables of Transcendental Functions 





Review 
Number Page 
57 123 
95 205-211 
95 205-211 
117 222 
56 121 
95 205-211 
95 205-211 
E 244 136 
53 81 
53 81 
70 133 
37 73 
57 123 
118 223 
37 73 
57 123 
37 73 
57 123 











TABLES FROM ANALYSIS 241 
Review 
Author /s Title Number Page 
Modified Bessel Functions of the Second Kind: Integral Order 
L. Fox A Short Table for Bessel Functions of Integer Order 37 73 
and Large Arguments 
W. FLUGcE Four-place Tables of Transcendental Functions 57 123 
Kelvin Functions of the First Kind 
W. FLicce Four-place Tables of Transcendental Functions 57 123 
Kelvin Functions of the Second Kind 
W. FLicce Four-place Tables of Transcendental Functions 57 123 
Kelvin Functions: Modulus and Phase 
E. JAHNKE, Tables of Functions E 243 42 
F. EMpE 
The Airy and Allied Integrals 
M. RoTHMAN Table of Bi(+-x) for x = 0.00(0.01)2.00and Bi(—x) 45 77 
for x = 0.00(0.01)10.00 
K. I. McKEnzieE, Table of Bi’(+-x) for x = 0(0.01)2 and Bi’(—x) for 69 133 
M. RoTHMAN x = 0(0.01)10 
Fresnel Integrals 
AxapemiiA Nauk SSSR = Tabliisy Integralov Frenelya [Tables of Fresnel 40 75 
Integrals] 
W. FLUcce Four-place Tables of Transcendental Functions 57 123 
Integrals involving Bessel Functions 
M. ROTHMAN Table of the integrals of Ai(+x) for 44 77 
x = 0.00(0.01)2.00 and x = 0.0(0.1)10.0 
M. RoTHMAN Table of the integrals of Bi(+x) for 46 78 
x = 0.00(0.01)2.00 and of Bi(—x) for 
x = 0.0(0.1)10.0 
TERENCE BUTLER, Tables of Definite Integrals involving Bessel Func- 50 79 
Kari POHLHAUSEN lions of the First Kind 
ASCHER OPLER, Tables for Predicting the Performance of Fixed Bed 51 79 
Nevin K. HIESTER Ion Exchange and Similar Mass Transfer Processes 
Cecit HASTINGS, JR., Approximations for Digital Computers 56 121 
JEANNE T. HAYWARD, 
James P. Wone, Jr. 
G. Eason, “On certain integrals of Lipschitz~-Hankel type 68 132 
B. NoBLE, involving products of Bessel functions” 
I. N. SNEDDON 
ALBERT D. WHEELON, A Table of Integrals Involving Bessel Functions 119 223 
Joun T. ROBACKER 
Complete Elliptic Integrals of the First and Second Kinds 
Cecit HAsTINGs, JR., Approximations for Digital Computers 56 121 
JEANNE T. HAYWARD, 
James P. Wona, JR. 











Author /s 


TABLES FROM ANALYSIS 


Title 


Incomplete Elliptic Integrals of the First and Second Kinds 


W. FLtcce 


S. RusHTON 


MiLTon ABRAMOWITZ, 
Partie RABINOWITZ 


Four-place Tables of Transcendental Functions 


Hypergeometric Functions 


“On the confluent hypergeometric function 
M(a, Y sy 


“Evaluation of Coulomb wave functions along the 
transition line” 


Debye Functions, Radiation Integrals, Fermi-Dirac Functions, etc. 


J. CLUNIE 


Ann T. NEEMS 


A. J. W. DutjveEstijn, 
J. BERGuuIs 


“On Bose-Einstein functions” 


“Graphs of the Compton energy-angle relationship 
and the Klein-Nishina formula from 10 Kev to 
500 Mev” 


The Computation and the Expansion of Some Triple 
Integrals Originating from the Theory of Cosmic Rays 


Miscellaneous Integrals Involving Circular and Exponential Functions 


L. BERGHUIS 


ADMIRALTY RESEARCH 
LABORATORY 


H. GELLMAN, 
JEAN TUCKER 
M. ROTHMAN 
M. RoTHMAN 
M. RoTHMAN 


M. RoTHMAN 


Cecit HAstTINcs, JR., 
JEANNE T. HAYWARD, 
James P. Wone, Jr. 


A Table of Some Integrals 
“Table of |Z| =|/0® secde™se%da|”” 


“Tables of the functions Do(x) and D;(x)” 


“Tables of the integrals and differential coefficients 
of Gi(+x) and Hi(—x)” ” 


Table of Gi(x) and its derivative for 
x = 0.0(0.1)25(1)75 


Table of Hi(—x) and its derivative for 
x = 0.0(0.1)25(1)75 


Tables of the integrals of Gi(x) and Hi(—x) for 
x = 0.0(0.1) — 20.0 


Approximations for Digital Computers 


Miscellaneous Generalizations of Common Functions 


Staff of the Computation 


Department of 
Mathematisch Centrum 


E. W. DiyKstTRA, 
A. Van WIJNGAARDEN 


HERBERT E. SALZER 
HERBERT E. SALZER 


Table of Polylogarithms, Report R24, Part I: 
Numerical Values 


Interpolation Coefficients and Formulas 
Table of Everett's Interpolation Coefficients 


Formulas for Calculating Fourier Coefficients 
Formulas for Inverse Osculatory Interpolation 





Revi 
Number Page 
57 123 
24 37 
25 38 
26 38 
30 40 
120 - 223-224 
28 39 
41 76 
42 76 
43 77 
47 78 
48 78 
49 78 
56 121 
29 40 
89 202-203 
93 204-205 
94 205 














TABLES FROM STATISTICS 


Author/s Title 
Coefficients for Integration and Summation using Ordinates 
H. E. SALzEr “A simple method for summing certain slowly 


convergent series”’ 


H. E. SALZER “Equally weighted quadrature formulas over semi- 
infinite and infinite intervals” 


Coefficients for Integration and Summation, using Differences 
Owen R. Mock Fifty Everett Integration Coefficients 


Curve-Fitting, etc. 


HERBERT E. SALZER Formulas for Calculating Fourier Coefficients 
E. S. PEARSON, Biometrika Tables for Statisticians 
H. O. HARTLEY 


Harmonic Analysis—Sines and Cosines 


L. W. PoLLak Eight-place Supplement to Harmonic Analysis and 
Synthesis Schedules for Three to One Hundred Equi- 
distant Values of Empiric Functions 


k 
OweEN R. Mock A Table of Sines and Cosines = 


TABLES FROM STATISTICS 


Tables with tie Normal Deviates as Arguments 


D. V. LInDLey, Cambridge Elementary Statistical Tables 
J. C. P. MILLER 
Cecit HASTINGS, JR., Approximations for Digital Computers 


JEANNE T. HAYWARD, 
James P. Won, Jr. 


W. FLicce Four-place Tables of Transcendental Functions 


E. S. PEARSON, Biometrika Tables for Statisticians 
H. O. HARTLEY 


Tables with the Normal Area as Argument, including 


percentage points and probits 
Cecit HASTINGS, JR., Approximations for Digital Computers 
JEANNE T. HAYWARD, 
James P. Won, Jr. 
E. S. PEARSON, Biometrika Tables for Statisticians 
H. O. HaRTLEY 
Rankings of Normal Samples or Order-statistics, Median 
H. RuBEN “On the moments of order statistics in samples 
from normal populations” 
E. S. PEARSON, Biometrika Tables for Statisticians 


H. O. HaRTLEY 


243 
Review 
Number Page 
91 203-204 
92 204 
90 203 
93 204-205 
95 205-211 
34 72 
79 196 
4 26 
56 121 
57 123 
95 205-211 
56 121 
95 205-211 
13 32 
95 205-211 








TABLES FROM STATISTICS 





Review 
Author/s Title Number Page 
Measures of Location and Measures of Normal Dispersion 
other than the Standard Deviation 
D. V. LInDLeEy, Cambridge Elementary Statistical Tables 4 26 
J. C. P. MILLER 
J. H. CapwELL “The statistical treatment of mean deviation” 10 31 
E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 
H. A. Davin, “The distribution of the ratio, in a single normal 105 216 
H. O. HARTLEY, sample, of range to standard deviation”’ 
E. S. PEARSON 
A. J. Duncan “The use of ranges in comparing variabilities” 107 217 
Max HALPERIN, “Tables of percentage points for the Studentized 111 219 
SAMUEL W. GREENHOUSE, maximum absolute deviate in normal samples” 
JEROME CORNFIELD, 
JuLta ZALOKAR 
Tests of Normality Criteria 
E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 
Miscellanea (Normal Distribution) 
A. C. CoHEN, Jr., “Tables of Pearson-Lee-Fisher functions of singly 7 29 
Joun WoopwarD truncated normal distributions” 
H. WEILER “A new type of control chart limits for means, 11 31 
ranges, and sequential runs”’ 
R. E. BECHHOFER “A single-sample multiple decision procedure for 16 34 
ranking means of normal populations with known 
variances” 
R. E. BECHHOFER, “A single-sample multiple decision procedure for 17 34 
Mitton SOBEL ranking variances of normal populations” 
C. W. DunneETT, “A bivariate generalization of Student’s ¢-distribu- 18 34-35 
MiILTon SOBEL tion with tables for certain special cases’”’ 
E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 
E. S. PaGE “Control charts for the mean of a normal popu- 103 215-216 
lation” 
S. B. CHowpHuRY “The most powerful unbiased critical regions and 114 220-221 
the shortest unbiased confidence intervals associ- 
ated with the distribution of the classical D*- 
statistic” 
H. J. GAWLIk A Table of a Function Related to the Error Function 117 222 
Area of Incomplete I’-function 
E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 
Ordinates of x*-Distribution and Poisson Frequencies 
PEARSON, Biometrika Tabies for Statisticians 95 205-211 














TABLES FROM STATISTICS 245 
Review 
Author/s Title Number Page 
Percentage Points of x*-Distribution, Percentage Points 
of Pcisson Series 
D. V. LinDLey, Cambridge Elementary Statistical Tables 4 26 
J. C. P. Miter 
E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 
T. Lewis 99.9 and 0.1% of the x?-Distribution”’ 108 217-218 
Laguerre Series, Gram-Charlier Series B, Non-Central x? Series, 
Fractional Degrees of Freedom 
S. H. ABpEL-Aty “Approximate formulae for the percentage points 109 218 
of the non-central x?-Distribution” 
Cumulative Binomial Series, Tail Area of I, (a, b) 
H. G. Romic 50-100 Binomial Tables 5 28 
LEs.iE E. Srmwon, Tables of the Cumulative Binomial Probabilities 36 72 
Frank E. Grusss 
Individual Frequencies of Binomial Distribution 
H. G. Romic 50-100 Binomial Tables 5 28 
RSMTC Table of Binomial Coefficients 31 70 
LEs.re E. Smon, Tables of the Cumulative Binomial Probabilities 36 72 
Frank E. GrusBBs 
E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HarRTLEY 
Percentage Points of J, (a, b), Confidence Belts for 
Binomial Distribution 
E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 
T. E. STERNE “Some remarks on confidence or fiducial limits” 104 216 
Tables of t-Distribution 
D. V. LinDLey, Cambridge Elementary Statistical Tables 4 26 
J. C. P. Mr_er 
J. M. SENGUPTA “Significance level of 2x*/ (2x)? based on Student’s 19 35 
distribution” 
E. S. PEARSON, Biometrika Tables for Statisticians 95 205-211 
H. O. HARTLEY 
MARTIN WEIBULL “The distribution of ¢- and F-statistics and of cor- 106 217 
relation and regression coefficients in stratified 
samples from normal populations with different 
means” 
Tables of F-Distribution and Multivariate Variance Tests 
26 


D. V. LInDLEy, Cambridge Elementary Statistical Tables 4 
J 








Author /s 


E. S. PEARSON, 
H. O. HARTLEY 


MARTIN WEIBULL 


H. G. Romic 


S. H. ABpeL-Aty 
J. C. Sprrz 
J. W. WHITFIELD 


Ws. R. THOMPSON 


E. S. PEARSON, 
H. O. HARTLEY 


MarTIN WEIBULL 


J. W. WuHITFIELD 


E. S. PEARSON, 
H. O. HarRTLEY 


H. A, Davin 


TABLES FROM STATISTICS 


Title 
Biometrika Tables for Statisticians 


“The distribution of t- and F-statistics and of cor- 
relation and regression coefficients in stratified 
samples from normal populations with different 
means” 


Hypergeometric Series 
50-100 Binomial Tables 


Miscellaneous Discrete Distributions 
“Tables of generalized k-statistics” 
“Matching in psychology” 


“The distribution of the difference in total rank 
value for two particular objects in m rankings of n 
objects” 


Tables of the Four Variable N- and Psi- Functions 


Distribution of Ordinary Correlations 


Biometrika Tables for Statisticians 


“The distribution of t- and F-statistics and of cor- 
relation and regression coefficients in stratified 
samples from normal populations with different 
means” 


Fisher’s r to z and z to r Transformations 
Cambridge Elementary Statistical Tables 


Biometrika Tables for Statisticians 


Rank Correlation and Associated Criteria 


“The distribution of the difference in total rank 
value for two particular objects in m rankings of n 
objects”’ 


Biometrika Tables for Statisticians 


Rankings, Gaps, etc., for Rectangular and other 


Non-Normal Distributions 
“The distribution of the range in certain non- 
normal populations” 
Order Statistics from any Continuous Parent 


“Universal bounds for mean range and extreme 
observation” 


“Statistical treatment of censored data. Part I. 
Fundamental formulae” 





Review 
Number Page 
95 205-211 
106 217 
5 28 
9 30 
20 35-36 
21 36 
22 36 
95 205-211 
106 217 
+ 26 
95 205-211 
21 36 
95 205-211 
116 222 
15 33 
96 211-212 














Author /s 
BENJAMIN EPSTEIN 


J. H. CaDWELL 


S. ROSENBAUM 
C. K. Tsao 


T. J. TERPSTRA 


F. G. Foster, 
A. STUART 


H. E. DANIELS 


YosHIHIKO HrRaco, 


HimeEnorr Morrmoura, 


Hisao WATANABE 
W. J. Drxon 


TABLES FROM STATISTICS 


Review 
Title Number 


“Tables for the distribution of the number of 99 
the exceedances” 


“The probability integral of rangeforsamplesfrom 112 
a symmetrical unimodal population” 
Non-Parametric Two Sample and Analysis of 
Variance Randomization Tests 
“Tables for a nonparametric test of location”’ 12 


“An extension of Massey’s distribution of the 100 
maximum deviation between two-sample cumula- 
tive step functions” 


Runs 
“Theory of the probability distribution of runs in 8 
a sequence of observations” 
Other Non-Parametric Tests 


“The exact probability distribution of the T sta- 14 
tistic for testing against trend and its normal 
approximation” 


“Distribution-free tests in time-series based on the 97 
breaking of records” 


“A distribution-free test for regression parameters” 98 
“Tables for three-sample test”’ 101 


“Power under normality of several nonparametric 102 
tests”’ 


Formulas and Tables of Pearson Type Frequency Distributions 


E. S. PEARSON, 
H. O. HARTLEY 


J. CLunre 


Cecit HASTINGS, JR., 
JEANNE T, HAYWARD, 
James P. Wonca, Jr. 


A. M. FREUDENTHAL, 


E. J. GuMBEL 
S. R. BROADBENT 


D. R. Cox 


W. L. STEVENS 


Biometrika Tables for Statisticians 95 


Tables of Miscellaneous Continuous Distributions 


“On Bose-Einstein functions” 26 
Approximations for Digital Computers 56 
“Minimum life in fatigue” 110 


“The quotient of a rectangular or triangularanda 113 
general variate”’ 


Distribution of Means from Non-Normal Parents 


“The mean and coefficient of variation of range in 115 
small samples from non-normal populations” 


Angular Transformations 
“Tables of the angular transformation” 6 


247 


Page 
213 


219 


32 
213-214 


30 


33 


213 


214 


214-215 


205-211 


38-39 
121-123 


218-219 


220 


221-222 


28 





TABLES FROM PHYSICAL SCIENCES 





Revi 

Author/s Title Number Page 

Random Numbers 
D. V. LINDLEY, Cambridge Elementary Statistical Tables 4 26-27 
J. C. P. MILLER 
Random Deviates from Mechanical Devices. 
Quasi-random Numbers 

K. D. TocHER “The application of automatic computers to 67 132 


sampling experiments” 


Random Normal Deviates 


K. D. TocHEeR “The application of automatic computers to 67 132 
sampling experiments” 


TABLES FROM PHYSICAL SCIENCES 


Note: Many tables whose calculation was directly instigated by research in 
physics are purely mathematical in nature and these are classified 
under Analysis or another appropriate heading. 


Mechanics of Solids 
James G. BERRY Tables of some Functions Related to the Legendre 27 39 
Functions P,-™(x) and Q,(x) when n is a complex 
number 
Fluid Mechanics 
AERONAUTICAL RESEARCH A Selection of Graphs for Use in Calculations of 52 80-81 
CouncIL Compressible Airflow 
K. G. TADMAN Tables of Flow Functions for Bodies of Revolutionin 53 81 
Circular Tunnels and Jets 
Cecit HASTINGS, JR., Approximations for Digital Computers 56 121 


JEANNE T. HAYWARD, 
James P. Wone, Jr. 


Exchanges, Reaction Kinetics, Reaction Rates, etc. 


ASCHER OPLER, Tables for Predicting the Performance of Fixed Bed 51 79 
Nevin K. HIeEsTER Ion Exchange and Similar Mass Transfer Processes 


Electricity and Magnetism 


CHESTER SNOW “Formulas for computing capacitance and in- 70 133 
ductance”’ 
Atomic Physics 
MriLTton ABRAMOWITZ, “Evaluation of Coulomb wave functions along the 25 38 
PHILIP RABINOWITZ transition line” 
Nuclear Physics and Scattering 
(Including Photon Scattering and Compton Effect) 
Ann T. NELMs “Graphs of the Compton energy-angle relationship 30 40 


and the Klein-Nishina formula from 10 Kev to 
500 Mev” 

















Author /s 


H. GELLMAN, 
JEAN TUCKER 


Cecit HAstincs, JR., 
JEANNE T. Haywarp, 
James P. Won, Jr. 


A. J. W. Dutyjvestiyn, 
J. Bercuuts 


J. CLunie 


AUTOMATIC COMPUTATION 


Title 
“Tables of the functions D(x) and D;(x)” 


Approximations for Digital Computers 


Cosmic Rays 
“The Computation and the Expansion of some 
triple integrals originating from the theory of 
cosmic rays” 


Statistical Mechanics—Classical and Quantum 


“On Bose-Einstein functions” 


AUTOMATIC COMPUTATION 


Conferences and Symposia on Computing Machines or Calculation 
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